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SYl^/IBOLS  USED 


A:  Plate  thickness  (Nondimens ional  length  unit) 

b:  A  half  of  the  layer  thickness 

N:  Layer  number  in  the  plate 

p:  Density 

t^:  Impact  time 

P^(x^,t):  Impact  stress 

*  Xj^(n),  X2(?):  Coordinate  variables 

*  t(T):  Time  variable 

T  :  Nondimens ional  time  unit  =  a//c~7p 
o  5o 

*  Stress  tensor  and  its  components 

*  u,,  u(U),  v(V):  Displacement  vector  and  its  components 

0 

*  c..-  ,  c..(C..):  Elastic  Moduli  (c--  =  c-- - 

ijk£'  ij  xj^  11  11  0^2 

A,  p:  Lame's  constants 

Infinitesimal  strain  tensor 

ij  1 

Laplace  transform  (in  t)  and  Fourier  transform  (in  n)  of  A 
P^(5):  Legendre  polynomial  of  g 

*  k(K) :  Wave  number 

*  (A)((i)):  Frequency 

0,  a,  g:  Phase  shifts  (wave  number  through  the  thickness) 

C^:  Dilatational  and  shear  wave  speed 
G*(cl))  =  G'(m)  4-  iG”(a)):  Complex  modulus  of  elastomer 
D:  Thickness  of  viscoelastic  layer 

h:  A  half  of  the  crack  length 

:k 

Quantities  in  (  )  are  nondimensional  quantities. 
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Preface  and  Summary 


This  report  is  the  last  of  a  series  on  the  response  of  composite 
plates  to  impact  forces.  The  motivation  for  these  studies  has  been  an 
attempt  to  understand  the  damage  to  aircraft  jet  engine  fan  blades  by 
foreign  object  impact  such  as  ice  balls,  stones,  and  birds.  In  addition, 
the  National  Aeronautics  and  Space  Administration,  sponsors  of  this 
research,  have  sought  to  develop  computer  codes  from  these  analyses  which 
will  aid  the  fan  blade  designer  in  locating  potential  failure  modes  and  posi¬ 
tions  and  thus  assist  in  optimizing  fan  blade  fabrication  to  create  greater 
impact  tolerance. 

The  basic  approach  of  the  principal  investigator  in  these  studies  has 
been  to  use  wave  propagation  techniques  to  model  the  early  response  of 
composite  plates  to  impact  type  forces  In  using  the  wave  method,  the 
plate  can  be  simplified  in  the  analyses  by  neglecting  reflections  from  edge 
boundaries  far  from  the  impact  point.  Thus, while  the  overall  geometry 
of  the  plate  is  no  longer  included  in  the  analysis,  more  sophisticated 
mathematical  models  near  the  point  of  impact  have  been  used. 

The  basic  model  for  the  composite  plate  studies  has  been  the  anisotropic  plate 
theory  as  extended  by  Mindlin  [  1  ]  to  account  for  wave  phenomena-  The 
plate  equations  were  used  as  an  approximation  of  the  exact  theory  of 
elasticity  because  they  lead  to  simpler  computational  schemes  for  evalua¬ 
ting  average  stresses  and  displacements  in  the  plate. 

Fourier  and  Laplace  transform  techniques  have  been  used  throughout 
these  studies  and  inversion  of  the  transforms  has  been  accomplished 
with  a  fast  Fourier  transform  algorithm.  This  algorithm  is  an  effective  compu¬ 
tational  tool  but  ^requires  careful  scaling  of  the  impact  problem  in  both  space 


and  time 


variables.  When  it  is  properly  used  it  can  lead  to  calculations  of  thousands  of 
stress  values  in  a  fraction  of  the  time  of  conventional  finite  element  schemes. 


In  summary,  the  use  of  plate  models  for  the  fan  blade  impact  has  avoided 

the  analytical  complexities  of  the  exact  theory  of  elasticity  as  well  as 
the  computational  difficulties  of  finite  element  methods. 

In  earlier  reports  both  central  and  edge  impact  of  an  anisotropic 
plate  were  studied,  [  2-4].  In  those  reports  only  wave  propagation  in 
the  plane  of  the  plate  was  investigated.  In  another  report  [5]  a 
multilayer  plate  model  was  developed  in  order  to  study  impact  induced 
wave  propagation  in  both  the  thickness  and  in  plane  directions.  In 
this  final  report  further  results  are  presented  from  the  multilayer 
model.  The  composite  plate  has  been  modeled  with  as  many  as  eight 
separate  layers.  Each  layer  may  itself  have  several  plys,  so  that 
effective  anisotropic  constants  must  be  used  for  each  layer  in  the  analysis. 

The  mathematical  model  exhibits  wave  propagation  in  both  the  thickness  and 
inplane  directions.  Impact  generated  waves  are  shown  to  lead  to  inter¬ 
laminar  shear  stresses  and  interlaminar  tensile  stresses  during  and  after 
impact . 

This  report  also  presents  an  analysis  of  an  impact  damping  mechanism. 

This  consists  of  thin  damping  layer  introduced  between  composite  layers  in  the 
mathematical  model.  The  resulting  response  due  to  impact  shows  that  considerable 
reduction  of  stress  can  be  achieved.  However  it  appears  that  this  stress 
reduction  is  linked  to  the  lower  elastic  moduli  of  the  damping  sublayers 

and  not  the  viscous  nature  of  the  sublayer. 

Fianlly  an  attempt  was  made  to  analyze  the  impact  of  a  plate  with  a 
crack.  While  the  problem  has  been  formulated,  no  progress  was  made  on 
obtaining  numerical  answers  to  the  crack  problem. 


IV 


I.  INTRODUCTION 


The  present  research  is  a  continuation  of  our  previous  work  on  the 
stress  wave  propagation  in  a  laminated  composite  [2-5].  It  has  been 
motivated  by  the  problem  of  the  impact  on  jet  engine  fan  blades  caused  by 
ingestions  of  foreign  materials,  such  as  birds  and  hailstones.  The 
successful  application  of  fiber-reinforced  composite  materials  depends  on 
the  ability  of  these  materials  to  withstand  forces  due  to  such  impact. 

The  simplest  approach  to  examine  the  dynamic  behavior  of  a  composite 
plate  is  based  upon  the  work  of  White  and  Angona  [6].  In  their  work, 
referred  to  as  the  effective  modulus  theory,  the  response  of  the  composite 
plate  to  waves  propagating  normal  to  the  laminate  is  predicted  by  a 
single  constant  wave  speed, regardless  of  the  internal  structure  of 
composites.  Even  though  this  theory  is  satisfactory  for  many  problems,  it 
fails  in  the  case  of  some  problems  when  the  wave  lengths  become  short .  To 
overcome  this  limitation.  Sun  and  et  al.  proposed  a  model  which  includes 
the  effects  of  internal  structure,  such  as  the  layer  thickness  [Tl*  Iti  their  work, 
referred  to  so  the  effective  stiffness  theory,  displacements  of  both  the  reinforo- 
ing  layer  and  the  matrix  layer  are  expressed  as  linear  expansions  about  the 

midplanes  of  the  layers  and  approximate  equations  of  motion  are  derived 
for  both  layers.  Then  these  approximate  equations  are  required  to  satisfy 
the  continuity  conditions  of  displacement  and  stress  components  on  every 
interface.  Using  this  model  the  propagation  of  harmonic  waves  has  been 
examined . 

More  recently  a  number  of  researchers  have  presented  models  for  multilayer 
plates  either  by  the  discrete-continuum  theory  or  the  continuum  mixture 
theory  [8-l4].  Many,  however  examined  only  the  frequency- 

wave  number  dispersion  relationship  and  stopped  short  of  the  transient 
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impact  problem  except  for  a  fev  experimental  or  n-americal  -works  using 
the  finite  element  method  which  sometimes  show  a  considerable  discrepancy 
from  the  experimental  results  • 

In  this  report  we  present  a  new  attempt  to  mathematically  model  the 
multilayer  plate  and  develop  a  method  by  which  we  can  examine  the 
transient  propagation  of  an  impact  wave  in  the  plate,  not  only  along 
the  longitudinal  direction  but  also  through  the  thickness  direction  of  the 
plate  as  well,  using  an  inexpensive  Fast  Fourier  Transform  algorithm 

[3,15]. 

The  composite  plate  under  consideration  for  the  first  part  of  the 
present  report  is  imagined  to  comprise  N  identical  elastic  layers.  And 
each  layer  is  made  of  a  number  of  unidirectional  plys  lying  alternately  at 
a  layup  angle  i(()  from  the  symmetry  axis^as  shorn  in  Fig,  1.  Then  the 
elastic  properties  of  the  plate  depend  on  the  layup  angle  (|).  A  key 
assumption  for  the  first  step  of  the  work  is  that  all  the  layers  are 
identical.  While  restricting  the  application,  this  assumption  allows  us 
to  formulate  the  problem  using  difference-differential  equations  due  to 
a  rather  simple  periodic  structure  of  the  plate.  This  technique  for 
periodic  structures  has  been  widely  used  in  the  study  of  electrical 
transmission  lines  [16]  and  in  the  vibration  of  multistory  buildings  [17]* 
By  means  of  an  approximate  plate  theory  of  Mindlin  [I8],  a  set  of  approx¬ 
imate  equations  of  motion  is  developed  for  a  typical  layer  using  the  inter¬ 
laminar  stresses  as  explicit  variables.  The  relative  motion  of  a  layer 
to  the  adjacent  layers  is  related  by  phase  shifts  which  represent  the 
solution  of  the  difference  parts  of  equations.  In  this  way  the  number  of 
the  layers  can  be  increased  without  increasing  the  size  of  matrix  in  the 
numerical  process  of  invert  to  satisfy  the  boundary  conditions. 
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It  is  also  well  understood  that  a  thin  viscoelastic  layer  present 
between  elastic  layers  can  reduce  the  elastic  impact  energy  significantly 
by  dissipating  the  strain  energy  into  heat  [19,20].  In  our  previous  work 
[5]  an  elastomer  layer  is  presented  between  a  composite  half  space  and  a 
protection  strip  on  the  edge  on  which  the  impact  is  applied.  Numerical 
results  of  the  work  showed  an  appreciable  reduction  in  the  normal  stress. 

As  an  extension  of  this  research  and  the  first  part  of  this  report 
we  now  examine  the  wave  propagation  in  a  composite  plate  made  of  two 
elastic  layers  and  an  elastomer  layer.  Generalization  of  this  problem  is 
straightforward  by  assuming  that  our  new  periodic  composite  layer  is  now 
made  of  an  elastic  sublayer  and  a  viscoelastic  sublayer  lying  alternately. 

W.e  can  now  develop  the  approximate  theory  which  includes  both  sub¬ 
layers.  For  the  second  part  of  the  present  research  we  will  examine  the 
simplest  case  of  this  kind,  i.e.,  an  impact  on  a  composite  plate 
consisting  of  two  elastic  layers  and  an  elastomer  layer  between  them. 

Another  possible  extension  of  the  multilayer  composite  plate  which 
can  be  found  in  frequent  practice  is  discussed  in  the  last  part  of  this 
report.  In  this  chapter  a  crack  is  introduced^  on  the  interface 
between  two  elastic  layers  which  represent  the  final  step  before  a 
failure  occurs  in  the  composite  either  by  spalling  or  by  excessive  shear 
stress.  Such  crack  problems  constitute  the  main  part  of  the  study  of 
fracture  mechanics.  A  serious  mathematical  difficulty  arises  even  in  the 
static  problems  because  of  the  mixed  boundary  conditions  along  the  crack 
direction.  The  difficulty  becomes  more  serious  in  the  case  of  dynamic  problems 
due  to  the  diffraction  of  waves  at  the  crack  tip  [21-24].  By 
employing  the  approximate  equations  of  motion  developed  in  the  first 
part,  the  transient  wave  problem  has  been  formulated  and  dual  integral 


equations  are  obtained  after  application  of  the  mixed  boundary  conditions. 

But  the  resulting  dual  integral  equations  are  not  easy  to  solve  and  are  under 
investigation  at  this  time. 

In  the  results  presented  in  this  report  only  a  line  impact  has  been 
examined.  This  has  simplified  the  calculations  and  saved  computer  time  in 
testing  the  model.  The  technique,  however,  can  be  extended  to  the  two- 
dimensional  or  central  impact  problem.  Since  the  impact  speed  is  very  high 
(-450  m/sec)  and  the  impact  time  is  short  (<  100  ysec) ,  the  impact  can  be 
in  the  range  of  the  elastic-plastic  impact  or  even  in  the  range  of  the 
hydraulic  impact.  But  the  initial  transmission  of  impact  energy  is  propa¬ 
gated  by  elastic  waves,  as  if  in  anunbounded  plate,  and  it  is  useful  to 
investigate  the  problem  by  means  of  the  linear  theory  of  anisotropic 
elasticity  in  an  infinite  composite  plate. 


II.  IMPACT  ON  MULTILAYER  ELASTIC  PLATE 


1.  Formulation 

Basic  Theory  of  Linear  Anisotropic  Elasticity 

Cauchy’s  equations  of  motion  in  cartesian  tensor  form  without  body 
forces  are  given  by 


ij  »i 


=  pu. 


a. .  =  a. . 

13  31 


(II-l) 


where  the  repeated  index  implies  summation  on  that  index*  A  comma  repre¬ 
sents  a  partial  differentiation  with  respect  to  the  index  following  the 
comma  and  a  superposed  dot  denotes  a  time  derivative, 
tensor  is  related  to  the  infinitesimal  strain  tensor  by 

-It 


a. . 
13 


or  a . 
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13 
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(II-2) 


The  condensed  elastic  moduli  c^^  has  the  following  form  for  orthotropic 
materials 
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Analysis  of  a  Layer 

For  a  layer  shown  in  Fig.  1  we  employ  the  approximate  plate  theory 
of  Mindlin  [l8  ]  and  the  displacement  field  u  is  expanded  in  terms  of 
the  Legendre  polynomials  as 


u(Xj^,X2,x^,t) 


I  u^^^Cxj^.x^.t)  •  p^(5) 
n=0 


(II-3) 


where  5  is  the  local  coordinate  along  ihe  thickness  direction,  normalized 
by  b,  a  half  layer  thick. 

Instead  of  solving  Eq.  (Il-l)  directly  we  solve  a  new  approximate 
equation  of  motion  which  is  obtained  through  a  variational  process  by 
integration  of  Eq.  (II-l)  over  the  thickness  C  (see  [I],-  [23]).  The  result  is 


b 


aja 


+  [P^CO.a^.l 


1 

C=-l 


2pb  ..(n)  .  j  =  1,2,3 
2n+l  j  •  a  =  1,3 

(II-4) 


where 


= 


?(n)  ;  dP(E)  „  ,, 


By  substituting  the  constitutive  relation  CII’"2)  ^or  the  displacement 


expansion  (II-3)  into  the  above  approximate  equations  of  motion,  we  can 

^  (0)  (0)  (0)  (1) 

find  governing  equations  of  motion  in  terms  of  u^  ’  ^2  ’  ^3  ’  ^1 

The  accuracy  of  this  approximate  theory  depends  on  how  many  terms  of  the 


~  T  - 


displacement  field  we  retain.  Since  the  complexity  in  formulation 
increases  rapidly  with  the  number  of  terms  included  we  keep  terms  only  up 
to  second  order.  Furthermore,  we  will  examine  harmonic  waves  propagating 


along  the  x^  and  directions  so  that  we  drop 


(n) 


terms  and  set 


T —  {  }  =  0.  Next  to  get  rid  of  the  undesired  coupling  with  higher  modes 

CfXq 

^  ••fZ)  -fZ) 

we  set  =  ii^  ■'  =  0.  Then  the  resulting  equations  are 


2bpSf> 


‘42-“22’  ■  ^'’‘’“2 


(0) 


(1) 


2b,  (1)  ^3  (2).  (0)^1  (1).  ^  ^  ^  ••(!) 

T^‘^66'^2,ll''b‘^66^.1^  ■  2^"l2'^l,l^b‘^22'^2  >  <^^22^22^  =  3^P^2 


^®22"'^22^  -  C22U2  ^)  -  0 


(II-5) 


where  the  sign  +  and  -  represent  the  stress  components  on  the  top  and 
bottom  surfaces  of  the  layer  under  examination,  i.e.,  at  C  “  il*  Here 
we  notice  that  the  first,  fourth,  and  last  equations  are  written  in  terms 
of  u^^^ ,  where  (n+i)  is  an  odd  integer  and  represents  the  thickness 
stretching  motion  (or  symmetric  motion) .  In  the  rest  of  the  equations  in 
which  (n+i)  is  an  even  integer  the  displacements  represent  the  flexual 
motion  (or  antisymmetric  motion).  Hence, this  process  has  decoupled  the 
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stretching  motion  from  bending  motion I  To  get  rid  of  the  2nd  order 

(2) 

modes  from  Eq.  (II-5)  we  solve  the  last  two  equations  for  U2  and 

(2) 

u^  ,  and  insert  them  into  the  remaining  equations.  Then  Eq.  (II-5)  can 
be  reduced  as  follows: 


2b(c^^u^  2bpu^ 


'^66^b“l,l'‘^2,il^  ^°22~^22^  2.hpn^  ^ 


(II-6) 


(1)  , 


^  .•(!) 
3tpui 


where  ^12/^22- 

Plate  Analysis 

In  view  of  the  Legendre  polynomial  expansion,  the  displacements  on 

the  both  sides  of  a  layer  can  be  written  as  uT  =  uf^^  i  uf^^  since  the 

111 

governing  equations  for  a  layer,  Eq.  (11-6),  only  include  terms  up  to 
the  first  order  of  expansion,  i.e.,  a  linear  expansion.  Remembering  that 
this  analysis  is  valid  for  any  arbitrary  layer  in  a  plate,  say  the 
nth  layer,  equation  (II-6)  can  be  immediately  written  as 


These  two  motions  are? of  course^coupled  through  the  boundary  conditions. 
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^12  1 

C.  ^  (u  +U  T  )  +  "  •'••  (v  -V  -  )  +'r  (t  -T  ^  ) 

11  n  n-1  ,11  b  n  n-1  ,i  b  n  n-1 


P  1  ^  “ 

n  n-1 


•  b  '12  ‘“n-^b-l  V  4  '22  ‘  VVl>*f 


p  (ii  -ii  . )  =  c.  -  (u  -u  , ) - (u  -u  1 )  ”  'T  (v  +v  . ) 

n  n-l'^  11  n  n-1^^1  ^2  n  n-1  b  66'  n  n-1  ^ 


+ -  (a  -a  ^ )  +  — (t  +t  - ) 

c„T  n  n-1  1  b  n  n-1 


(11-7) 


P(Vr,+V  ,) 

n  n— 1 


=  c 


66 


{^(u  -u  - )  +(v  +v  - )  }  +  T ) 

b  n  n-1^1  n  n-1  ^11  bn  n-1 


where  a  and  T  are  used  to  represent  022  ^  ^ 

denote  and  U25  respectively.  These  equations  are  the  approximate 

equations  of  motion  of  a  layer  written  in  the  form  of  a  difference-differential 
equation.  For  a  plate  made  of  N  layer,  the  above  equations  contain 
4(N+1)  unknowns  and  offer  4N  equations.  Since 

the  additional  four  conditions  are  supplied  by  boundary  conditions  on 
the  top  and  bottom  surfaces,  solutions  of  these  equations  can  be  found. 

In  Eq.  (II-7)  we  notice  some,  important  points.  The  first  point  is 
that  the  logitudinal  coordinate  and  time  variable  are  continuous 

variables  while  the  thickness  coordinate  X2  is  now  discrete.  This 

enables  us  to  use  integral  transforms  in  x^  and  time  variables  so  that 
we  can  arrive  at  pure  difference  equations  after  integral  transforms. 

The  second  point  concerns  the  continuity  conditions  of  stress  and  dis¬ 
placement.  We  note  that  u,  v,  -^22  ’  ^12  continuous  across 

the  layer  boundary  and  these  conditions  are  identically  satisfied  by 
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Eq.  (II-7)-  But  the  normal  stress  tangential  to  the  layer  hoimdary  is  not  nec 

essarily  continuous  and  Eq.  (II-7)  allows  such  a  possibility.  One  can 
retain  higher  order  terms  in  the  displacement  expansion  given  by  Eq. 

(II-3)  to  give  more  accurate  results.  This  can  be  achieved  more  easily 
by  using  Eq.  (II-7)  and  increasing  the  number  of  layers  in  a  plate  under 
consideration.  This  process  does  not  give  any  additional  difficulties 
except  a  little  more  computer  time. 


2.  Dispersion  Relationships  of  Harmonic  Waves 
Harmonic  Waves 

Before  we  examine  the  transient  propagation  of  stress  wave  due  to  an 
impact  we  first  investigate  dispersion  relations'  of  harmonic  waves  in  a 
composite  plate  governed  by  approximate  equations  of  motion  (II-7) .  For 
harmonic  waves  propagating  along  the  axis  we  assume 


{u  ,v  ,a  ,T  } 
n  n  n  n 


{U  ,V  ,E  ,T  }  e 
n  n  n  n 


i(kx^-a)t) 


(II“8) 


Substituting  this  into  the  approximate  equations  of  motion  (II-7)  we  obtain 


-3c  iK(U  +U  )+(u  -3c„)  (V  -V  )+iKb(T  -T  )+3b(i:  +Z  O  =  0 

1/  n  n~l  22  n  n-1  n  n-1  n  n-1 


3c  iK(u  -U  )+(w^-c, (V  +V  )+b(r  -S  ,)  =  0 

00  n  n-1  DO  n  n~l  n  n-1 


(II-9) 


+  ~  iKb(Z  -Z  )  = 

Ci2  n  n-1 


0 
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for  n  =  1,  2  ...N.  Here  we  set 

K  =  bk  =  A  =  2bN 

-2  ,  2  2  2,  A  2 

0)  =  pb  0)  =  poj  (■^) 


and  A  is  the  total  thickness  of  the  plate.  For  a  plate  consisting  of 

N  layers,  the  boundary  conditions  require  traction  free  surfaces,  namely, 

=  0.  When  these  conditions  are  applied  to  equation 

(II-9)  we  obtain  4N  equations  in  terms  of  4N  unknowns  (U^,V^; 

U  ,V  ,T  ,1  with  n  =  1,...N-1;  U_-,V__)  .  By  setting  the  coefficient  matrix 
n  n  n  n  N  N 

to  be  singular,  required  dispersion  relationships  can  be  obtained. 


One-layer  Plate 


The  dispersion  relationship  for  a  plate  made  of  a  single  layer  can 


be  found  by  setting  N  -  1  in  equation  (II-9)  with  Z  =T  =Z-=T-=0. 

o  o  ±  1 

The  resulting  equations  are  now  written  in  matrix  form  as  follows: 


(II-IO) 


Then  by  setting  the  determinant  of  the  coefficient  matrix  zero  we  obtain 
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21-2  -2  2 
CiiK  -  — (o)  -3c22)(‘^  )  =  0 

2  1,-2  ^  2  -  .  ,-2  2,  ^ 

C6gK  -  3((.  -c^^K  -3Cgg)(co  -c^gK  )  =  0  . 


(II-ll) 


Here  we  notice  that  the  first  relationship  corresponds  to  the  state  of 

deformation  of  U-  =  U  and  V-  =  -V  ,  which  represents  the  thickness 
1  o  1  o 

extension  of  the  plate  (or  the  symmetric  mode),  and  the  second  describes 

the  f lexual  deformation  (or  antisymmetric  mode) .  The  exact  theory  of 

plates  gives  an  infinite  number  of  dispersion  relationships,  but  because  this 

model  only  has  two  inertia  points  (namely  n  -  0,  l) ,  each  of  them  having 

two  components  of  displacement,  we  only  have  the  first  four  relationships. 

Dispersion  relationships  and  corresponding  phase  velocities  for  an 

isotropic  plate  with  Poisson’s  ratio  1/4  (namely  X  =  p)  are  given  in 

Fig.  2a  and  2b  up  to  the  range  where  the  wave  length  becomes  equal  to  the 

plate  thickness.  Solid  lines  represent  the  symmetric  modes  and  dotted 

lines  the  antisymmetric  modes.  As  predicted  by  Mindlin  and  Medick  the 

optical  branch  of  the  S3rmmetric  mode  approaches  the  dilatation  wave  [l8]. 

The  acoustic  branch  of  the  antis3nnmetric  mode  starts  from  the  bending 

motion  and  approaches  the  shear  wave  when  the  wave  number  k  becomes 
* 

larger  and  larger  .  Similar  relationships  for  an  anisotropic  plate  made 
of  55%  graphite  fiber-epoxy  matrix  with  a  layup  angle  of  45°  are  shown  in 
Fig.  3a  and  3b. 


*  See  section  5  for  discussion  about  the  large  wave  number  limit. 
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Two-layer  Plate 

In  this  case  we  obtain  eight  equations  by  putting  n  =  1  and  2  in 
equation  (II-9) .  Boundary  conditions  require  T  =  E  =  T„  =  E-  =  0.  By 

o  o  ^ 

following  the  same  procedure  we  find  the  dispersion  relations  as 

r,-2  2^,-2  ,  ^  ,  2  2^,-2  -  2  ,  .  ‘^66*^12  2, 

{(0.  -c^^K  )(m  -3c22)-3c^2^  }(m  -c^^k  66"^“^  ^ 

+3((j)^-c^^k^){  (u^-CggK^)  (a)^-c^^K^-3cgg)-3Cggic^} 

..-2  2.  .-2  -  2  _  ,  -  2  2, 

+3(m^-CggK^){(m^-Cj^^K^)(m^-3c22)-3c^2'^^^  =  0 

Again  the  first  equation  represents  the  syimetric  mode  and  is  shown  as 
solid  lines  in  Fig.  4  and  5.  The  second  equation  is  plotted  with  dotted 
lines  representing  the  antisyininetric  mode. 

As  expected  we  have  six  relationships  since  the  this  two-layer  model  is  equiva¬ 
lent  to  a  three-mass  system  with  two  degrees  of  freedom  for  each  mass.  I^Jhen  the 
wave  number  kA  increases  the  following  are  observed:  for  the 

symmetric  mode  the  upper  optical  branch  approaches  the  dilatation  wave^vhereas 

* 

for  the  antisymmetric  mode  the  lower  optical  branch  approaches  the  shear  wave  . 


=  0 
(11-12) 


* 


See  section  5  ^'or  discussions  about  the  large  wave  number  limit. 
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N-Layer  Plate 

-2 

In  general,  we  can  obtain  a  2(N+1)  order  polynomial  of  o)  by  expand¬ 
ing  a  (4n)  X  (ta)  determinant  and  finding  2(N+1)  dispersion  relationships. 

But^  unfortunately^  this  process  involves  considerably  complicated  algebra 
and  it  may  be  necessary  to  develop  a  computer  technique  to  find  roots  of 
an  equation  in  determinant  form  (not  in  polynomial  form) • 

A  difference  equation  approach  can  be  used  to  solve  the  N  set  of 
four  simultaneous  first  order  difference  equatior^s  given  by  Eq.  (II-9)  . 

This  proceedure  is  neat  and  can  be  generalized  for  any  number  of  layers 
as  discussed  in  the  next  section; but  the  last  step  of  this  approach,  where 
a  long  polynomial  is  to  be  solved  again,  is  not  any  simpler  than  the  previous 


direct  method. 
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3.  Impact  on  an  Elastic  Composite  Plate 

Normalization  and  Integral  Transforms  of  Governing  Equations 

The  governing  equations  given  by  (II-7)  are  first  nondimensionalized 
as  follows: 


tu  =  {u  /A,U  /A,Xt/A} 

n  n  n  n  i 


'  ''=lj''=66-V''66'V‘=66^ 


T  =  t/T 

o 

where  A  is  the  total  thickness  of  the  plate  and  is  the  time  required  for 

the  quasi-shear  wave  to  travel  the  impact  radius .  Next  we  apply  a  Laplace 
transform  in  x  and  a  Fourier  transform  in  r|9  i.e.. 


g(s)  =  f  g(t)e 


i(k)  =  -i-  /  g(n)e^^’^dn 

-00 


Then  the  resulting  equations  are 


-(fs^  +  -:^  k^)(U  +U  J-CT.ik(V  -V  ,)+(T-T  .) 

2N  n  n-1  12  n  n-1  n  n-1 


=  0 


C  ik(U  +U  )-C^+2NC  )(V  -V  0  -T  i)  =  0 

Iz  n  n-1  3  22  n-1  n-1  n  n-1  6N  n  n-1 


(11-13) 


-C,.ik(5  -5  )-(fs^+-^  k^)(V  +V  )+(Z  -!  )  =  0 

DO  n  n-1  2N  n  n-1  n  n-1 


2c  k 

-(^+  +2NC-,)(U  -U  ,)+C  ik(V+V  (S-S  t)+(T-T  0=0 

3  6N  66  n  n-1  66  n  n-1  ^^^22  ^  ^ 
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where  the  normalization  factor  f  is  given  as 


f  = 


1  A  p 


66  o 


bAp 

c  ,T^ 
66  o 


Solution  of  Difference  Equations 

Since  the  simultaneous  difference  equations  given  are  linear 
and  all  the  coefficients  are  constants  the  solution  [26]  has  to  he 


A  A 


{U  ,V  ,T  =  {A,B,C,D}e 

n  n  n  n 


2i0ri 


(11-14) 


where  the  phase  shift  9  is  complex, in  general, and  propagation  through  the 

thickness  direction  in  the  plate  is  characterized  hy  6-  Namely,.  6  is  the 

wave  number  in  the  thickness  direction.  By  substituting  solution  (11-14) 

into  the  difference  equation  (11-13)  we  obtain  a  set  of  four  simultaneous 

homogeneous  equations  through  which  the  relationships  among  the  constants  A,B,C, 

and  D  have  to  be  determined.  If  we  set  the  resulting  coefficient  matrix  of 
A,B,C,  and  D  to  be  singular  we  obtain  the  following  equation  for  phase  shift  0: 


4  2  ^11^  2 

cos^9(fs^+-^^)(fs^-h- 


2N  ^ 


4  f  s 
+  sin  0(-  ® 


^12^  fs 


2  2  2 
+  cos  9  sin  9[ (f s  • 


— 1- 2NC 

^^^66  3  ^^^^22  ^6r  (^21  ' 


2  2 
11  2  12 


(11-15) 


=  a, 


4 

cos  9 


2^ 

a^cos  0 


+  a^ 


0  . 


This  equation  implies  that  for  a  given  wave  number  k  along  and  a 

frequency  s  s  represents  the  frequency  for  the  case  of  harraonic  waves), 
an  infinite  value  of  wave  numbers  exists  for  propagation  through  the 

thickness  direction, but  only  four  of  them  are  sufficient  to  give  all  linearly 
independent  solutions  of  the  form  of  Eq.  (11-14).  If  we  denote  the  solution 
of  the  phase  shift  equation  as 


cos  6  = 


(11-16) 


the  complete  general  solutions  of  difference  equation  (11-13)  are 


(11-17) 

Next, by  substituting  the  above  solutions  into  the  original  difference 
equations  (11-13)  we  find  the  relationships  among  A^,  C^,  and  . 

The  results  are 


Mi>  Hl>  <l>  cr!l> 
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X^(B) 


X^(3) 


X2(e) 


X3(e)  E3 


• cos2n3+i 


XjCB)  Ej 


XjCb) 


sin2n3 


(11-18) 


Y^(a)  E^ 


Y^Ca)  E3 


^2(0)  E3 


Y2(a)  E^ 


• cos2na+i 


‘  sin2na 


YgCa)  E^ 


Y3(a)  E3 


where  E-  =  D-  +  D^,  =  D-  -  D^,  E^  =  +  C, ,  E,  =  -  C,  and 
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X^(6) 


A.(3) 


.  /ON  2  ,  ^11,  2^  ,,  2  ,  ^66,  2^  3„ 

A(3)  =  (fs  +  +~2^k  )cos  3 


C  2  C 

+  {  (fs^  (^|— +-^^^+2nC,  ,)-C,,C,  _k^-c3ck^}sin^3'COs3 

2n  3  on  66  66  12  66 


A^ce) 


C  C 

=  ik  sin^8cos8{;'  (fs^  +  -;^^^)-(C  „+C,,)} 
6nc22  2n  12  66 


(11-19) 


A2(3) 


.  .  3„r^66^12,  2  ,fs^  ,  *^11,  2.,  ^  2„  .  2.^11,  2, 

=  1  Sin  3 {■■— ”— --k  --(—z — +2nC-^)}-cos  3sin3(fs  +  k  ) 
onc^o  3  on  66  2n 


2  C 

A3  (3)  =  ksin^6{(i|-  +  -g^2^2nCg^)C 


7  2 

^12^. 


3  -6n-  —66^-12-6^1^66^ 

,  .  3ln3cos^(^(.3^  .^)  .  (13^  .^kVds^  .^k^)C^,> 
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and 


Y^(a) 


A^(a) 

A(a) 


2  .  ^66,  2,  .fs^  , 


A(a)  =  i  cos^a*sina{Cggk^  (C  *  ^^^3  *"  6n 


^  k^  2 

+  2nCgg  +  6^)}  +  i  sin\(^+2nC22) 


^  ^  6nC22  "3  6n  ^66^ 


T  ^  N  3  ..  2^  66  ,  2, 

Ai(ot)  =  cos  a(fs  +"^  k  ) 


^  ,  2  ,  .k  2^12  2^^66,2.  .k^^ 

+  sxn  a.cosa{-(^)  — (fs  +-^k  ) +^Cgg+(— +  2nC22) > 


A2(a) 


k-cos  asina(C 


(11-20) 


.  ,  .  3  rl  ,fs^  .  ^  ^  ^12^66, 

+  ksin  aiT— (—7~+ — Z — +  2nC,,)-(-7— )  — - - ) 

6n  3  6n  66  6n  C22 


A^Ca) 


ik  sin^acosnC^^  2nc,6) 


2  C 

+  c  cos^a{-Cj^2(^®^  ^ 
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Equations  (II~18)  with  (11-19,20)  and  the  phase  shifts  a  and  3 
given  by  (11-16)  constitute  the  final  form  of  the  general  solutions  of  the 
difference  equations  (11-13) . 

In  Eqs.  (11-19,20)  we  notice  that  when  k  0  we  have  ~  ^^(S)  = 

¥2(01)  =  ^^(a)  =  0,  Namely  the  propagation  of  the  normal  stress  (with  phase 
shift  3)  and  the  propagation  of  the  shear  stress  (with  phase  shift  a)  are 
completely  decoupled.  This  occurs  when  the  waves  are  propagating  only  through  the 

thickness  direction  [27 ] . 

Impact  Boundary  Condition 

Boundary  conditions  for  an  impact  can  be  described  by  any  two  conditions 

among  u  ,  v  ,  a  ,  and  t  and  another  two  conditions  from  u„,  v„,  a._,  - 

®  o  o  o’  o  N  N  N  and 

For  our  present  problem  we  have  chosen  a  line  impact  by  a  normal  stress 

along  the  axis  (Figure  i. 

p 

00  =  — ^  (1-cos  '^^)  (1+cos^)  :  |x|  j<  a  and  0  £  t  £  to 

o 

=  0:  |x[>a  or  t>0  or  t>t 

'  ‘  o 

Hence, the  boundary  conditions  for  the  present  impact  problem  lead  to  "^he 


following  equation 
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0 

1 

1 

y  0 

0 

iX2(3)sin2  3N 

,  X2(B)cos23N  , 

cos2c!N 

cos23N 

,  isin23N  , 

iY2(a)sin2o(N 

.  0 

_ J 

0 

,  Y^(a) 

• 

"3 

q 

,  isin2ctN 

^3 

0 

,  Y2(a)cos2cxN 

"4 

0 

(11-22) 


where  q  is  the  integral  transform  of  the  impact  function  (II~21)  . 
Solving  the  above  equations  for  E^s  we  can  have 


V  =  {1+X2(B)Y^(oi)  }sin2aN.sin2gN  +  (.$)Y2  (a)  (.cos2aN*cos26N-l) 

*  X2(3)Y2(ot)  (cos2aN*cos23N-l)  +  sin2aN.sin26N 

^*2  “  i{cos23N.sin2aN-X2(e)Y2(a)coa2aN'Sin2aN}  (11-23) 

^3  ~  iX2(3){X^(e)Y2(a)sin23N-cos2aN-cos26N-sin2aN} 

*^°®23N‘cos2aN— 1}  • 


Substituting  the  E^s  into  the  general  solution  (ll-l8)  we  can  find 

the  complete  solutions  which  satisfy  the  impact  boundary  conditions  given 
by  (11-21).  In  other  words,  for  given  values  of  k  and  s  we  first  find 
the  phase  shift  a  and  3  from  (11-15,16)  and  with  these  we  can  find 
solutions  in  integral  transform  from  equations  (11-18,23)  x<rhich  are  the 

•A.  JS, 

final  solutions  under  impact.  After  U  ,  V  ,  f  and  Z  are  calculated 

n  n  n  n 


Allowing  the  determinant  of  the  coefficient  matrix  to  vanish  leads 
to  dispersion  relations  of  an  N-layer  plate,  namely  i?(a,3  )  =  0.  Then, 
a  and  3  are  obtained  from  ( 11-15, l6)  which  gives  the  complete 
dispersion  relationships. 
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with  a  given  impact  function  q,  they  can  be  inverted  easily  by  means  of  the 
fast  Fourier  transform  routine  [3,20]  to  give  the  complete  displacement  and 
the  stress  fields  after  impact. 

Tangential  Normal  Stress 

As  discussed  following  Eq.  (11-7), the  tangential  normal  stress  does  not 
appear  explicitly  in  the  approximate  equations  of  motion.  Therefore,  this 
component  of  the  stress  has  to  be  calculated  from  the  constitutive  equation. 
Namely, 


^11 


^  12  .  . 

■=11%,  1  -IT  (V''o> 


O,,  =  c,,u  ,  + 


'12 


11 


11  n,l  2b 


(v-v  ,)  l<n<N 
n  n-1  —  — 


or  after  normalization  and  integral  transform  they  are 


o 


(11-24) 


^11 


-ikC,,U  +  N.C,_(V  -V  .) 
11  n  12  n  n-1 


1  <  n  <  N 


Then  once  the  displacement  field  is  computed  the  tangential  normal  stress 
can  be  obtained  from  the  above  equation  and  inverted. 


4.  Some  Numerical  Results 


The  analysis  discussed  in  the  previous  section  includes  the  tiensient 

propagation  in  all  directions  but  suitable  choices  of  impact  time,  impact 
radius,  sizes  of  time  and  distance  steps  are  essential  to  make  good  use 
of  the  fast  Fourier  transform.  For  example,  if  we  take  a  large  time 
increment  with  a  relatively  thin  plate  propagation  through  the  thick¬ 
ness  will  not  be  seen.  For  this  matter  we  have  examined  several  different 
cases. 

Case  1;  Longitudinal  propagation 

Propagation  of  impact  generated  waves  along  the  longitudinal  direction 
is  examined  for  an  isotropic  plate  (steel  plate:  A  =  y  =  1.2x10^  psi) 
employing  a  two-layer  model.  For  these  calculations  we  used  an  impact 
time  t^  =  10  ysec,  plate  thickness  A  =  1  cm,  and  impact  radius  a  =  4  cm. 
Some  of  the  results  at  a  few  different  time  sequences  are  shown  in  Fig. 

6  a-f. 

In  these  figures  we  can  see  two  distinct  states  of  propagation  and 
corresponding  wave  fronts:  one  for  horizontal  displacements  (u)  and 
longitudinal  normal  stresses  another  for  vertical  displacements 

(v)  and  shear  stresses  (t).  In  other  words,  the  initial  signals  of  the 
horizontal  displacements  and  longitudinal  normal  stresses  propagate 
through  the  plate  with  longitudinal  wave  speed  at  amplitudes  that  are 
relatively  small.  But  the  major  parts  of  their  signals  are  due  to  a 
bending  wave  propagating  with  shear  stresses  and  vertical  displacements 
with  a  lower  velocity.  When  the  group  velocities  of  these  waves  are 
calculated  from  the  numerical  results,  they  are  about  5  mm/ysec  and 
3  mm/ysec,  respectively,  while  the  phase  velocities  of  the  unbounded 
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medium  of  this  material  are  =  / (A4’2|i) /p  =  5.61  mm/ysec  and 

C  =  /y/p  =3.25  mm/ysec. 
s 

Case  2:  Propagation  Through  Thickness 

To  examine  the  propagation  through  the  thickness  it  is  necessary  to  have 
a  sufficient  number  of  layers  in  a  plate.  It  is  also  essential  to 
make  the  time  step  relatively  small  compared  to  the  layer  thickness.  To 
do  this  we  increase  the  thickness  of  the  plate  and  the  number  of  layers  and 
reduce  the  impact  time. 

In  Figs. 7  and  8  propagation  of  the  transverse  normal  stress  in  the 
same  plate  (A  =  4  cm,  t^  =  2ysec,  a  =  40  cm;  4-layer  model)  is  shown 
at  different  time  sequences.  As  seen  in  Fig.  lithe  transverse  normal 
stress  is  initially  compressive  due  to  the  impact  and  a  compression  wave 
propagates  through  the  thickness.  But  later  it  becomes  a  tension  wave 
after  reflection  from  the  free  surface  and  the  tension  wave  propagates 
back  to  the  impact  surface.  In  Fig.  8  we  see  the  change  of  the  transverse 
normal  stress  and  the  interlaminar  shear  stress  with  time  for  the  same 
impact  conditions  as  in  Fig.  7. 

Similar  results  are  also  shown  for  the  case  of  an  anisotropic  plate 
in  Fig.  9  (55%  graphite  fibers-epoxy  matrix,  layup  angle  =  15°;  A  =  1  cm, 
t^  =  2ysec,  a  =  2  cm;  8-layer  model).  Here  we  again  notice  a  clear  delay 
in  time  for  waves  to  travel  from  one  layer  to  the  next  one.  Another  important 
point  is  that  the  shape  of  the  impact  stress  is  relatively  well  preserved 
during  the  initial  stage  of  propagation  but  changes  immediately  afterwards. 

The  distortion  of  the  shape  becomes  more  serious  with  further  propagation 
due  to  reflection,  thus,  showing  the  highly  dispersive  nature  of  the  harmonic 

waves  in  the  approximate  plate  theory. 


-  25  - 


When  the  group  velocities  are  calculated  from  these  results,  ve  find  6.32 
min/ysec  for  the  dilatation  -wave  and  3.33  mm/ysec  for  the  shear  wave  in  the 

case  of  the  isotropic  plate  and  2.5  mm/ysec  for  the  quasi-dilatation  wave  of 
the  anisotropic  plate.  Their  expected  values  are, respectively ,5.61, 

3.25, and  2.36  mm/ysec.  In  other  words,  waves  going  through  the  thickness 
are  traveling  faster  than  expected. 

Case  3.  Wave  Surfaces 

In  the  previous  two  cases  we  examined  the  transient  waves  propagating 
dominantly  along  either  the  x^  axis  or  through  the  thickness  direction 
by  suitable  choices  of  the  plate  geometry  and  impact  condition, 
now  examine  the  combined  effect,  simultaneous  propagation  in  both  direc¬ 
tions.  This  effect  is  shown  in  Fig.  10  (isotropic  plate;  A  -  4  cm, 
t^  =  4ysec,  a  =  4  cm;  4-layer  model)  where  the  transverse  normal  stress 
generated  from  the  line  source  due  to  impact  not  only  spreads  out  in  all 
directions  but  also  reflects  from  the  free  surface. 

When  the  plate  is  anisotropic,  the  situation  is  more  complex  in  the 
sense  that  waves  are  neither  dilatation  nor  shear  but  they  are  coupled 
together  (now  called  quasi-dilational  or  quasi-shear  waves) .  Due  to  the 
coupling, phase  velocities  of  the  anisotropic  wave  vary  from  one  direction 
to  another,  resulting  in  complicated  shapes  for  the  velocity  surfaces  and 
wave  fronts  [2].  For  an  ansiotropic  plate  (made  of  55%  graphite  fiber- 
epoxy  matrix  with  layup  angle  45^)  the  velocity  surfaces  and  the  wave 
surfaces  are  shoxm  in  Fig.  11.  The  stress  state  at  10  ysec  after  the 
impact  on  the  same  plate  (A  =  4  cm,  t^  =  4ysec,  a  =  2  cm;  8-layer  model) 
with  the  corresponding  wave  fronts  are  shown  in  Fig.  12a.  In  the 
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propagation  of  the  quasi-longitudinal  wave  we  notice  that  the  longi¬ 
tudinal  propagation  is  well  bounded  by  the  quasi-dilatational  wave  surface 
but  the  transverse  propagation  is  not.  The  shear  wave  is  not  bounded  by 
the  quasi-shear  wave  front  in  either  direction. 

This  interesting  phenomenon  of  higher  propagation  speeds  through  the 
thickness  is  related  to  the  dispersion  relationship  at  short  wave  length 
limits;  it  is  discussed  in  the  next  section. 


5.  Correction  Factor  and  Conclusion 


Correction  Factor 

According  to  the  present  model  of  a  multilayer  plate,  one  of  the 
antisymmetric  modes  of  the  dispersion  relation^ips  approaches  the  shear  speed 
when  the  wave  length  becomes  shorter  and  shorter,  as  mentioned  in  Section 
2.  It  is  well  understood  that  such  a  limit  is  incorrect,  i.e.,  in  the 
limit  of  short  wave  length  there  should  be  a  Rayleigh  wave  instead  of  a 
shear  wave.  Such  a  discrepancy  can  be  eliminated  by  introduction  of 
proper  correction  factors, as  shown  by  Mindlin  and  Medick  [l8  ] .  Correction 
factors  can  be  found  by  examining  either  the  large  wave  number  limit  or 
the  cut-off  frequencies  of  both  the  exact  theory  and  the  present  approximate 
theory.  Since  these  two  ways  lead  us  to  the  same  results  we  will 
find  the  factors  by  matching  the  cut-off  frequencies  of  the  two  theories. 

The  cut-off  frequencies  of  the  exact  theory  for  an  isotropic  plate 
can  be  obtained  from  the  well-known  Rayleigh-Lam’s  equation.  The 
lowest  cut-off  frequencies  are  ~  i/7x+2y ) /p  for  the  symmetric  mode  and 
~  V  y  /p  for  the  ant is3mme trie  mode.  The  corresponding  cut-off  frequencies 

2  r _  2  _ * 

of  our  approximate  theory  are  -g  ^^22^^  A*  ^ ^^66^^  ‘  notice 

2  2 

that  replacing  ^7  ^22*^  ^66  ^66^  makes  the  two 

theories  have  the  same  two  lowest  cut-off  frequencies.  Furthermore  the 
shear  wave  observed  in  the  short  wave  length  limit  of  the  present  approx¬ 
imation  becomes  a  wave  with  a  speed  of  — ^  i.e.,  the  Rayleigh  wave. 

/12 

Another  important  consequence  of  the  correction  factor  is  to  reduce 
propagation  speeds  through  the  thickness,  which  are  related  to  ^^22^^ 

* 

The  lowest  two  cut-off  frequencies  are  found  from  Eq.  (II-ll)  and  they 
are  independent  of  the  layer  number  in  the  plate  under  investigation. 
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and  }/c,Jp ,  with  a  factor  of  Propagation  of  the  maximum  value 

DO 

of  the  interlaminar  normal  stress  through  the  thickness  is  examined  with 
and  without  correction  factors  and  the  results  are  shown  in  Fig.  13. 

Without  the  correction  factor  the  propagation  speed  in  a  composite  plate 
is  roughly  about  2.60  mm/psec  obtained  from  the  numerical  results  used  in  Fig.  12. 

When  the  same  plate  is  subjected  to  identical  impact  conditions  this 
reduces  to  about  2.41  mm/ysec  with  the  correction  factor.  Comparing  this 
with  the  group  velocity  in  an  unbounded  composite  space  C~  2.36  mm/ysec) 
the  agreement  of  the  present  approximate  theory  is  remarkable.  Similar 
results  are  also  observed  in  "^he  case  of  shear  and.  (^uasi” shear  waves .  When 
these  correction  factors  are  introduced  in  previous  cases,  shown  in  Figs* 

8,  9^ and  129  all  the  signals  propagating  through  the  thickness  are  now 
well  bounded  within  the  corresponding  wave  fronts, as  shown  in  Fig.  12b 
and  from  this  we  can  notice  the  importance  of  the  correction  factors. 

Discussion  and  Computation  Time 

It  is  interesting  to  compare  the  computation  time  of  this  model  with 

some  other  methods,  such  as  the  finite  element  method  or  the  finite  difference 

method.  In  the  case  of  an  8-layer  anisotropic  plate  model,  from  which  Figs.  9 

and  12  are  produced,  we  have 

9  steps  along  the  thickness:  8-layer  model; 

32  step  along  the  x^  direction:  64  points  are  used  in  pratice  but 

only  half  of  them  are  useful  because 

of  the  symmetry  of  the  problem, 

32  steps  in  time; 

2  displacement  components  at  each  point . 

Therefore  the  total  number  of  the  unknowns,  which  are  the  basic  unknowns 
either  in  case  of  the  finite  difference  or  finite  element  methods,  is 


18,432.  After  these  primary  variables  are  calculated,  27,648  secondary 
variables  (three  stress  components  at  each  points)  have  to  be  calculated 
^Qcording  to  our  present  model  all  these  processes  require  only 
200  K  of  computer  space  without  using  magnetic  tapes  or  any  kind  of 
additional  storage  space  and  only  1  minute  6  seconds  for  CPU  time  in  the  IBM 
370-168  model  including  compiling,  linkage  editing,  I/O  and  execution. 

Conclusion 

The  present  theory  is  a  generalization  of  Kindlin' s  approximate  plate 
theory  applied  to  a  multilayer  plate  under  an  Impact.  By  combined  use  of  the 
finite  difference  technique  in  the  thickness  direction  and  the  fast  Fourier 
transform  in  the  plane  of  plate  and  time,  this  model  can  be  very  useful 
for  the  study  of  wave  propagation  in  a  composite  plate  under  impact  forces. 
However,  reasonable  attention  in  usage  of  the  fast  Fourier  transform  is 
required  to  avoid  spurious  data.  From  the  limited  numerical  data  obtained 
from  this  model  it  appears  that  the  anisotropy  in  the  plate  will  lead  to  a 
considerable  interlaminar  shear  which  might  lead  to  ply  bonding  failures. 

The  model  also  shows  that  for  short  enough  impact  times,  an  interlaminar 
tension  can  develop  as  one  would  expect ,  which  might  also  account  for 


interlaminar  ply  failure. 


111.  IMPACT  OF  A  COMPOSITE  PLATE  WITH  M  INTERLAMINAR  DAMPING  LAYER 


1,  Description  of  Problem 
Geometry  of  Plate 

As  an  extension  of  the  multilayer  plate  discussed  in  Chapter  II 

ve  nov  examine  the  impact  and  the  consequent  stress 

wave  propagation  in  a  composite  plate  with  viscoelastic  damping  layers. 
Possible  models  for  damping  mechanisms  in  plates  are  shown  in  Fig.  14. 

We  will  formulate  a  model  made  of  an  alternating 

number  of  elastic  and  viscoelastic  layers,  as  showTi  in  Fig.  l4-c. 

As  long  as  the  layer  structure  of  the  plate  is  periodic,  the  main  part  of 
the  analysis  in  Chapter  II  for  an  elastic  plate  is  valid  with  additional 
equations  for  viscoelastic  layers. 

Viscoelastic  Property  of  Elastomer 

The  mechanical  properties  of  an  elastomer  are  usually  expressed  in 
terms  of  a  complex  modulus  depending  on  the  frequency,  i.e., 

G*(a))  =  G'(aj)  +  iG”(a3)  .  (III-l) 

With  this  the  constitutive  equation  is  written  as 

a  (u)  =  G  (o))e  ((d)  (III-2) 

in  the  frequency  space  where  a..(a))  and  e..((D)  are  respectively  the 

13  ij 

'ft 

Fourier  transforms  of  a.,  and  e.,  in  time  [291. 

ij 

* 

For  (III-2)  the  Fourier  transform  is  defined  as 

00 

f((D)  =  f  f(t)e~^“^dt  . 

/2ir 
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The  constitutive  relation  (III-2)  with  the  complex  modulus  (III-l) 

implies  the  following  constitutive  i^uation  in  a  time  space: 

t 

a(x,t)  =  f  G(t-T)e(x,T)dT  (III-3) 

— .00 

where  the  relaxation  function  G(t)  is  related  with  the  complex  modulus 
* 

G  (m)  as 

00 

G*((o)  =  G(t)e"^“^dt  .  (III-4) 

o 

Therefore^ when  the  complex  modulus  G  (to)  is  obtained  by  experiments, 
usually  by  means  of  harmonic  excitation  of  strain,  the  relaxation  function 
G(t)  can  be  found  by  inversion  of  equation  (111^4). 

The  viscoelastic  property  of  the  elastomer  under  consideration  has 
been  extensively  investigated  [19])  and  its 

complex  modulus  is  given  in  Fig,  15.  This  complex  modulus  can  be  reasonably 
well  described  by  a  three  parameter  equation  as 

G'Cm)  =  a  -  -y—  .  (III-5) 

0)  +c 

These  three  parameters  are  obtained  from  another  set  of  parameters:  the 

maximum  values  of  G*(a))  when  the  maximum  value  of  G”(a))/G’ (m) 

and  the  at  which  G"(aj)/G' (w)  becomes  the  maximum.  Therefore  >if 

we  characterize  the  complex  modulus  by  proper  choice  of  G’Cm),  o)  and 

o 

the  maximum  of  G”(a3^)/G' (o)^)  ,  the  relaxation  functions  are  completely 


described. 


2.  Formulation 


Elastic  Layer 

In  Fig.  16  a  typical  viscoelastic  layer  (nth)  is  shown  between  two 
adjacent  elastic  layers  (nth  and  (n+l)th)  with  appropriate  discretiza¬ 
tion.  The  approximate  equations  of  motion  for  the  nth  elastic  layer 
given  by  Eq.  (II-7)  are  still  valid.  But  remembering  the  new  discretizing 
notation  in  Fig.  16  we  now  have  to  replace  (  )  and  (  )  -  by  (  ) 

and  (  ’  ^respectively .  The  results  are 


p(u  ) 

n  n-1 


/  - ,  +  V  ,  12  /  -  +  >.  il/“  +  \ 

c-  T  (u  -Ki  T )  T  T  +-^(v  -V  -T  T  ) 

11  n  u-1  ,11  b  n  n-1  ,1  b  n  n-1 


p 

n  n 


3c  3c 

^^12,-.+  ,  '"''22,  +  X  .  3,  , 

b  n  n-X  ,1  n  n-1  b  n  n-1  n  n-1  ,1 


p(u  -u^  !> 


3c  3c 

-  /  -  +  s  66,  -  +  ^  ^^66,  4-  , 

c-  -  (u  -u  1  )  T  T - ;^(u  -u  -  ) - r— (v  4v  t  i 

11  n  n-1  ,11  ,2  n  n-1  b  n  n-1  ,1 


(III-6) 


,11,-  +  .  _L  3  ,  +  V 

c  ’ 


P(v"+v^  .) 

n  n-1 


1^  T+(v +v'^  t)  n}+^(a  -a"*”  -) 

66  b  n  n-1  ,1  n  n-1  ,11  b  n  n-1 


Viscoelastic  Layer 

Since  the  thickness  of  the  elastomer  is  thin  compared  with  the  elastic 

layer,  we  can  assume  that  the  stress  field  is  uniform  through  the  thickness 

-  +  —  + 

of  the  elastomer.  In  other  words, we  have  a  =  a  ==  cr  and  x  =  x  =  t 

n  n  n  n  n  n 

for  (III-6) .  Therefore, the  following  compatibility  conditions  for  the 
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elastomer  can  be  obtained  immediately: 


=12(n)  -  4  2D<V"n’ 

(III-7) 

=22(t.)  - 

where  D  is  the  thickness  of  the  elastomer. 

We  further  assume  that  the  dissipation  is  mostly  due  to  shear  motion, 
i.e.,  that  the  normal  component  of  the  continuous  traction  vector  is 
transmitted  through  the  viscoelastic  layer  purely  elastically.  There¬ 
fore,  by  combining  (III-7)  with  (III-3)  we  find 


o 

n 


E.  +  -- 

— (v  -V  ) 
D  n  n 


T 

n 


/  G(t-T)  {j  j 

2  n  n  D  n  n 


(III-8) 


These  two  equations  and  four  more  from  Eq.  (III-6)  are  the  complete 
equations  needed  to  solve  the  impact  on  a  composite  plate  with  elastomer. 
For  a  plate  made  of  N  elastic  layers  and  (N^l)  viscoelastic  layers 
Eq.  (III-6)  provides  4N  equations  and  (III-8)  gives  2(N-1)  equations. 


Since  the  total  number  of  the  unknowns  are  now  6N+2  (u  ,  v  ,  a  ,  t  ; 

o  o  o  o 


Ui,  v^,  u^,  v^,  T^; 


’  '^N-l’  ^N-1’  ^-1’  \-l’  S-1’  '^N-1’ 


v^,  T^)  we  can  solve  this  system  of  equations  with  four  additional 

conditions  supplied  by  the  suitable  boundary  conditions. 


Here  we  notice  that  the  governing  equations  are  now  a  set  of  six 
dif f erence-integro-partial  differential  equations.  These  equations  can  be 
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reduced  to  difference  equations  after  appropriate  integral  transforms  and 
the  resulting  difference  equations  can  be  handled  rather  simply,  as  in  the 
previous  chapter. 

3.  Numerical  Results  and  Discussion 
Impact  on  Plate 

For  the  report  ve  examine  the  impact  on  a  plate  consisting 

of  two  elastic  layers  and  a  viscoelastic  layer, as  shown  in  Fig.  17 ^with 
an  impact  function 

X-j  2 

a  =  p  {l-( — )  }  sin  —  ;  Ix-  I  <  a  and  0  <  t  <  t 

o  oa  t’'l‘  o 

o 

(III-9) 

=  0  ;  Ixj^l  >  a  or  t  >  t^  ,  or  t  <  0 

with  all  other  stress  components  vanishing  on  both  surfaces  of  the  plate. 

Now  by  putting  n  =  1  and  2  into  Eq.  (III-6)  we  have  eight  equations  and 
two  more  equations  are  obtained  from  Eq.  (III~8) .  We  again  normalize 
these  equations  and  take  the  integral  transform,  as  in  Chapter  II.  The 


resulting  equations  are : 
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-(fs^+-4^bk^)(U*+U  )-C,„ik(v"-V  )+T_  =  0 

A  1  o  Iz  i  o  1 

3C  A 

3Cj2ik(0-^C^)-(fs2+-2i-)(f-4__)+3!j-ifc|f^  =  -3f^ 


-(fs^+^,.k^+^3C.,)(U,-U  )+3C,,ik(v7+V  ) --:^  YikE,+3T, 
A  11  b  66  1  o  66  1  o  0^2  All 


kd 


=22  *  “ 


A  I  A 


-(fs^+-|ibk^)(U2+U^)-C^2i’^(V^i^"^l  “  ° 


(III-IO) 


-(fs2+|£^lk2  3.|3C5j)  (02-O^)+3Cj,jik(f2+^)  |lk!i+3f^ 


=  0 


eA(V+-Vi) 


G(s){|(u^-u~)  --yCv^+v^)} 


where  G(s)  is  the  Laplace  transform  of  the  relaxation  function  G(t) 


with  respect  to  t  =  t/T^  and  we  have  used  the  boundary  conditions 
Tq  =  T2  =  cr2  =  0  .  From  the  above  10  equations  we  can  find  10  unknowns 
(U^,  V^,  U^,  V^,  U^,  V^,  E^,  T^;  U2?  ^2^  with  given  impact  function 
and  once  these  are  calculated  the  displacement  and  the  stress  fields  can 
be  computed  by  inversions  of  the  integral  transforms  by  means  of  the  FFT 


alogorithm. 
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Numerical  Results 

For  the  present  computation  we  have  used  the  Young’s  modulus 

12 

A  42  41*10 

E  =  ,7*10  psi  and  the  shear  modulus  G’(aj)  =  .817*19 - ^ — r — ^  for 

3*10+0) 

the  elastomer  where  ca  is  given  in  hertz.  The  G'(w)  in  this  case  implies 
that  G'(“>)  =  .817*10^  psi  and  max(G"(a))/G' (to))  =  3.3  at  “  ^00  Hz. 

The  propagation  of  stress  wave  in  this  case  is  quite  similar  to  that 
of  the  composite  plate  without  an  elastomer  layer  except  the  peak  values 
of  the  interlaminar  stress.  Values  of  the  peak  stress  with  different 
thickness  of  the  elastomer  layer  are  plotted  with  those  of  the  purely  elas¬ 
tic  plate  in  Fig.  l8.  As  we  can  see  in  this  figure  the  interlaminar  shear 
stress  has  increased  by  a  small  amount  while  a  reduction  of  the  normal 
stress  is  considerable  when  the  elastomer  layer  becomes  thicker  and  thicker. 
From  this  result  it  is  obvious  that  the  reduction  of  the  normal  component 
of  stress  can  be  achieved  by  introducing  such  a  soft  and  energy-dissipating 
elastomer  layer. 

Discussion 

In  addition  to  the  simple  reduction  of  the  normal  stress  it  is  also 

observed  that  the  amount  of  reduction  increases  with  the  value  of  G”(a3)/G’ (oj) 

and  the  location  of  o)  at  which  G"(a))/G’ (o))  becomes  the  maximum  value. 

o 

In  other  words,  we  can  make  the  dissipation  effect  more  serious  by  choosing 

an  elastomer  whose  G”(a))/G’ (o))  becomes  maximum  at  o)  around  which  the 

o 

most  of  the  impact  energy  is  carried  out. 

It  is  also  believed  that  a  further  dissipation  effect  will  be  possible 
if  we  make  the  transmission  of  the  normal  stress  viscoelastic  across  the 
elastomer  layer,  which  we  have  assumed  is  elastic  for  this  report. 


IV.  IMPACT  ON  A  PLATE  WITH  A  CRACK 


1.  Introduction 

When  the  impact  stress  is  low,  the  impact  is  elastic  and  the 
stresses  in  the  plate  can  be  described  by  elastic  wave  propagation.  When 
the  stress  is  increased  beyond  a  certain  limit  then  the  impact  damage 
occurs.  Elastic-plastic  impact  is  complicated  for  two  reasons,  namely, 
unloading  and  loading  must  be  treated  differently,  and  the  strain  rate  effect 
[30]  must  be  included.  If  the  impact  stress  is  increased  further  to  a 
certain  level  where  the  induced  stress  is  higher  than  the  strength  of  a 
target  material  then  penetration  begins  to  occur.  In  this  limit  the 
target  material  sometimes  behaves  as  a  fluid  and  such  a  state  of  impact  is ' 
known  as  a  hydraulic  impact  [31].  Another  failure  mode  is  the  occurrence  of 
interlaminar  cracks. 

Investigation  of  the  stress  state  in  solids  with  cracks  falls  in  the 
category  of  so-called  fracture  mechanics  and  has  been  under  an  extensive 
scrutiny  since  the  famous  enunciation  by  Griffith  [32].  Presence  of  cracks 
inside  a  material  usually  leads  us  to  a  mixed  boundary  value  problem  and 
only  a  limited  class  of  problems  can  be  solved  [33,34].  In  the  case  of 
dynamic  loading  the  problem  becomes  more  difficult  due  to  the  scattering 
of  the  stress  wave  by  the  crack  [21-24] .  In  this  report  we  will 
formulate  the  problem  of  a  plate  with  a  crack  which  is  subject  to  a 
dynamic  loading. 

Our  original  goal  was  to  study  the  effect  of  interlaminar  cracks  in 
composite  plates  in  response  to  impact  loads.  Debugging  problems  in 
other  parts  of  this  report, however, used  valuable  time  originally  set  aside 
for  this  problem.  The  following  section  is  an  attempt  to  illustrate  the 
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use  of  the  Mindlin  plate  theory  for  the  study  of  interlaminar  cracks  and 
to  point  out  the  mathematical  difficulties  that  must  be  overcome  in  solving 


the  problem. 


2,  Formulation 


Description  of  Problem 

The  plate  under  consideration  has  a  crack  on  the  midplane  running 
from  x^=  “h  to  4*h  as  shown  in  Fig.  19.  Stress  can  be  applied  either 
on  the  surface  of  the  plate  or  on  the  crack  surface.  In  the  former 
case  the  crack  surfaces  can  be  in  contact  and  the  boundary  conditions 
become  more  complex  due  to  the  partial  continuity  of  stresses  and  dis¬ 
placements  during  the  contact.  For  the  present  report  to  illustrate  the 
mathematical  difficulties  we  assume  that  the  crack  surface  is  subject  to 
a  known  compressive  impact. 

Governing  Equation  and  Boundary  Conditions 

We  can  formulate  this  crack  problem  by  assuming  that  the  lower  and 
the  upper  half  plates  are  made  of  a  number  of  layers  but  for 
simplicity  we  consider  the  plate  to  consist  of  two  identical  layers  and  the 
crack  to  be  present  on  the  interface  of  these  two  layers.  Following  the 
notation  shown  in  Fig.  19  we  have  the  governing  equations  identical  to 
Eq.  (III-6)  with  n  =  1  and  2.  The  boundary  condition  requires  that 
both  plate  surfaces  remain  traction  free.  The  crack  surface  is  subject 
to  a  prescribed  impact  condition  while  the  displacement  and  stress  are 
continuous  along  the  layer  boundary  outside  the  crack.  Namely, we  have 


a  "  T 


=  02  =  -^2  =  ° 


01  = 

+  — 

T,  =  T,  =  0 


x|  <  h 


u,  =  u. 


cr,  =  a. 


,  <  =  V- 


X  >  h 


(IV-1) 


-  39 


-  1+0  - 


Due  to  the  twofold  syiranetry  of  the  problem  we  now  have 

+  —  +  —  +  “ 

Vq  =  "V^,  =  0  and  we  can  set  ”  ^1’  ”^1  “  ~  ^1’ 

Thus^  the  eight  equations  obtained  from  Eq.  (III-6)  are  now  reduced  to 


p(u^+u^)  =  Cii('^l+“o\ll'^~^V^o\l 


p(VVo) 


p(u  -u  ) 
1  o 


^^1 2  ^^22  3o 

— •'T 


(IV-2) 


"ii^V%^  ,11  -—^^1^0^ 


P(^l-^o)  “ 


and  the  boundary  condition  is  now 


a  =  -P^(x^,t)  |xj  <  h 


Vl=  0 


X  >  h 


along  X2  =  0 


(IV-3) 


Dual  Integral  Equation 

We  now  normalize  the  governing  equation  (IV~2)  and  take  the  integral 
transform.  Then  we  have 


.  fc  2  .  3A^  . 

-3C^2^^  ,  (fs  +“g^22^’ 


U,+U 
1  o 


V  -V 
1  o 


0  ,  (fs^|c^^k2+fCg,),  -3Cg,ik 

_  L  _  L  J 


(IV-4) 


and  these  can  be  solved  for  U  ,  U- ,  V  ^  and 

o  1  o 

the  mixed  boundary  conditions  are  given  by  a 


in  terms  of  Z.  Since 
and  v^  we  solve  as 


=  K(s,k)E 


(IV-5) 


with 


A  =  Det 


-3C^2^k 


,  (fs  +-^22^ 


(IV-6) 


B  =  Det 
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Next  we  take  the  inverse  transform  of  Z  and  V^,  and  apply  the  mixed 
boundary  condition  given  in  Eq.  (IV-3) .  Since  the  boundary  conditions  are 
for  all  times  t  >  0  we  only  take  the  inverse  Fourier  transform  to  apply 
the  boundary  conditions,  i.e.. 


2(n,s)  =  —  /  I  e"^^’^dk 

i/JtT  -oo 


V^(n,s)  =  —  I  K(s,k)  Z  e'^^’^dk 

t/2tt  -00 


(IV-7) 


Application  of  the  boundary  condition  given  by  Eq.  (IV-3)  results  in  the 
following  integral  equation: 


-  1+2  - 


-P  (n,s)  =  —  /  l  e  :  Ini  <  h/A 

^  /Jrr  -a> 


0  =  —  /  K(s,k)  S  e'^^'^dk  :  |nl  >  h/A 

-<» 


for  an  unknown  function  E. 


(IV-8) 


3*  Discussion 

The  integral  equations  of  the  type  given  in  Eq.  (IV-8)  are  known  as 
dual  integral  equations, each  of  which  has  its  own  region  of  application 
and  occur  in  mixed  boundary  value  problems  [35],  There  are  a  number  of 
ways  to  solve  this  type  of  integral  equations,  such  as  by  reduction  to  a 
single  Fredholm  integral  equation  or  by  using  the  Wiener-Hopf  technique  [36] . 
Finding  the  solution  depends  on  the  kernel  and  in  general  it  is  rather 
difficult  to  do  except  for  some  special  cases  such  as  for  Bessel  kernels  or 

trigonometric  kernels , 

Once  the  unknoTxm  function  Z  is  determined  the  other  variables 
(Uo,Ui,Vo,Vi)  can  be  computed  by  solving  the  algebraic  equation  (IV-4) 
and  the  complete  displacement  can  be  found  by  inversions  of  'the  integral 
transform. 

The  problem  formulated  in  this  chapter  is  the  simplest  impact 
problem  in  that  the  contact  of  the  crack  surface  does  not  occur  and  that 
it  has u  twofold  symmetry.  But  it  is  expected  that  the  critical  response 
of  the  plate,  particularly  the  stress  field  near  the  crack,  can  be  a  guide“ 
line  for  a  more  complex  problem. 


V.  COUCLUSIOK  AUD  RECOMMEMDED  RESEARCH 


The  present  theory  is  a  generalization  of  Mindlin's  approximate 
theory  of  plate  applied  to  a  multilayer  plate  Tinder  impact.  By  combined 
use  of  the  finite  difference  technique  in  the  thickness  direction 
and  integral  transforms  this  model  has  been  shown  to  be  very  effective 
for  wave  propagation  analyses . 

This  model  is  extended  to  examine  the  effects  of  an  elastomer  layer 
between  elastic  layers  of  the  plate.  The  reduction  of  interlaminar  normal 
stress  is  signficant  due  to  the  damping  layer  but  further  investigation 
seems  necessary  to  determine  the  nature  of  the  reduction. 

The  presence  of  a  crack  in  the  plate  has  been  formulated.  The  re¬ 
sulting  equations  are  given  by  dual  integral  equations  which,  as  in  many 
cases,  are  rather  difficult  to  solve. 

The  basic  idea  of  the  periodic  structure  of  the  multilayer  plate, 
where  the  governing  equations  are  derived  for  each  layer  and  given  by  a 
set  of  difference- differential  equations,  may  be  useful  to  handle  dif¬ 
ferent  types  of  problems, such  as  heat  conduction  and  thermoelastic 
problems  in  composite  plates. 
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Dispersion  Relationship  and  Phase  Velocity  of  Composite  Plate 
One-Layer  Model  (55%  Graphite  Fiber  -  Epoxy  Matrix;  Layup  Ang 


etric  Mode 


Lgure  5.  Dispersion  Relationship  and  Phase  Velocity  of  Composite  Plate 
Two-Layer  Model  (55%  Graphite  Fiber  Epoxy  Matrix,  Layup  Ang 


Longitudinal  Propagation  in  Isotropic  Plate;  2-layer  Model 


Continued 


Continued 


Impact  surface 


Figure  8.  Transverse  Propagation  of  Normal  and  Shear  Stress 
(Same  as  in  Figure  7) 


Figure  9  Transverse  propagation  of  normal  stress  in 

composite  plate;  8-layer  Model  (55%  Graphite 
Fiber  -  Epoxy  Matrix,  ±15°  Layup;  A  =  1  cm, 
t  =2  usee,  a  =  2cm) 
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X2(mm)  Surface  (1  fJLSQC  after  impact) 


Figure  11.  Velocity  Surface  and  Wave  Surface  of 
Composite  Plate  (55%  Graphite  Fiber- 
Epoxy  Matrix,  Layup  Angle  45®) 
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Figure  12a  .Wave  front  10  fxsec  after  impact  (without  correction  factor) 
(557o  Graphite  Fiber-Epoxy  Matrix,  ±45°  Layup;  A  =  4  cm, 

8- layer  Model;  t  =  4  ysec,  a  =  2  cm) 
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Figure  12b.Wave  front  10  fisec  after  impact  (with  correction  factor) 
(55%  Graphite  Fiber-Epoxy  Matrix,  ±45°  Layup;  A  =  4  cm, 
8- layer  Model;  t=4ysec,  a=2  cm) 
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FIGURE  15b.  THE  COMPLEX  YOUNG’S  \;ODULUS  OF  LR3-604  MEASURED 
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Figure  19.  Composite  Plate 


APPENDIX  A  FLOW  CHAET 


In  this  flow  chart  and  program,  U(I,J),  V(I,J),  TAU(I,J) 
SIGMA  (I,J)  and  SIGMA1(I,J)  represent  U,  V,  f  and  E  in 
Eq.  (11-17,18)  and  integral  transform  of 
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Block  A  Inversion 


(i)  Data  xx(I,J)  in  integral  transformed  space  are  generated 
for  a  half  of  the  Inverse  transform  range:  I  =  1  ~  NT, 

J  =  1  -r  NX/2 

(ii)  Generate  full  data  over  J  =  1 ~  NX  by  FLIP 

Symmetric  flip;  V(I,J),  SIGMA(I,J),  SIGHAl(I,J),QO(i,j) 

Antisymmetric  flip:  U(I,J),  TAU(I,J) 

(iii)  Invert  them  for  displacement  and  stress  fields  by  FOURT 

(iv)  Take  care  of  the  coordinate  shifts  and  multiplication 
factors  in  FOURT  by  FACT 

(v)  Print  out  by  MAP 
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APPENDIX  B  SAMPLE  COMPUTER  DECK 


Data  Card  4 


6 • E—  Ub 


Data  Card  3 


Data  Card  2  I! 


15  0.c:456E 


0.400UE  06 


0.1170E  Or 


355EE  Ob 


Data  Card  1  |j 


'■GD.SVS'IH  BD  * 


/  PROGRAM 

/  P'/  EKEC  FORTGCLG 


PROGRAM  IN  SOURCE  DECK  (see  Appendix  C) 


I  i‘l  I TS  REG  I  DH=2  0  r  K  5  L I  HeS=5K  ?  CLASSIC 


JOB  CARD 


1  2  3  4  5  6  7  8  9  tfl'll  >2  ^3  M  13  IS  n  II  «  M  71  22  23  24  25  25  27  28  29  30  31  32  33  34  35  35  37  38  39  4  0  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  53  55  60  SI  62  63  64  65  £6  67  68  69  70  71  72  73  74  75  78  H  78  79  BQ 

n  1 1 1 11 1 1 1  n  11 1 1 11 1 1 1 11 1 1 1  n  1 1 1 11 11  n  11 11 1 1  n  n  11 1 11 11  n  n  11  n  n  11  n  m  n  1 11 11  n 
22222222 2 2222222222222?2222222222222222222222222222222222?2 2 22222222222222222222 
33333333333333333333333333333333333333333333333333333333333333333333333333333333 
44444444444444444444444444444444444444444444444444444444444444444444444444444444 
55555555555555555555555555555555555555555555555555555555555555555555555555555555 
68S668G6666866666G6SS66BS666B666666G6B666686668686fi66E886GG666B6665B66BB66B6666B 


8  8  8  8  S  S  8  8  8  8  8  8  8  a  8  8  8  8  8  8  8  8  8  8  8  8  8  5  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  S  3  8  8  ?  8  a  8  8  8  8  8  8  B  8  8  b  8  8  8  8  fi  6  8  6  c 

9999999999®9BS99999^9999939999999939899999939995S9S99999399999S9SS98S99SS9999999 
I  2  3  4  5  6  7  1  '213:1  17  16  13  Li  25  2-^2623  30  3-:  3^  43  M  .5  46  47  «  4^3 -.2  53  64  55  16  57  S3  SB  60  61  62  53  6- 65  66  67  66  69  70  71  T:  73  74  75  76  77  71  16  SC 


APPENDIX  C  LISTING  OF  PROGRAM  MD  SAMPLE  OUTPUT 


RELEASE  2.3 


MAIN 


DATE  =  77139  21/11/39 


C  « j|;  *  5|:  5|t  J|!  4: sjt  *  sit  5|s  ^  3|:  5ls  sic  sit «  Jj:  si' *  ***  **5!= ’Is ’!=***=!«***«***=!' * 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PROGRAM  CALCULATES  THE  TRANSIENT  PROPAGATION  OF  STRESS  WAVE 
IN  A  LAMINATED  COMPOSITE  PLATE  DUE  TO  A  NORMAL  IMPACT. 


THE  PRESENT  PROGRAM  IS  A  PART  OF  THE  RESEARCH  PROJECT  OF 
PROFESSOR  FRANCIS  C.  MOON 

DEPARTMENT  OF  THEORETICAL  AND  APPLIED  MECHANICS 
CORNELL  UNIVERSITY 

AND  SUPPORTED  BY  NASA-LEWIS  RESEARCH  CENTER. 


r  -T*  or*  or 


r  or  *r  ■rr  or  or  -O'*  O' 


I#  .Jt.  JU  .JU  «Cr  %)U  WU  «(■.  «L>  »U  .A#  .JU  Wt.  WU  Wk.  ^  ^  <4*  jJc  siS  SSs  5*i 

l■»3^^.^»•.^«.^■nroro'ororo'^oro'or^o'oro''r'f'oror^o'ororo'o*oro*o*■T•o'•'ro«'r»■o“0'o*ororo'0'^'* 


FOLLOWING  USER'S  GUIDES  ARE  PROVIDED  BY  DR.  B.S.KIM  AND  ALL  THE  EQUATION 
NUMERS  CORRESPOND  EQUATIONS  IN  THE  ACCOMPANYING  TECHNICAL  REPORT 


1.  DATA  TO  BE  SUPPLIED  BY  4  DATA  CARDS  (IN  THE  ORDER  OF  READING) 


lANGLE:  FIBER  LAYUP  ANGLE  IN  COMPOSITES 

C11tC12..  :  ELASTIC  MODULI  OF  COMPOSITE  LAYERdN  PSD 

NLAYER:  NUMBER  OF  THE  LAYERS  IN  THE  GIVEN  PLATE 
DfDELTA;  THICKNESS  OF  COMPOSITE  PLATEdN  CM) 

RHO:  MASS  DENSITY  OF  COMPOSITE  LAYERdN  GR/CM*=:'3) 

NA:  RADIUS  OF  THE  IMPACT  AS  A  MULTIPLE  OF  THE  PLATE  THICKNESS 
TAUO:  IMPACT  TIMEdN  SECOND) 

NX:  NUMBER  OF  INEGRATION  STEPS  IN  SPACE 
NT:  NUMBER  OF  INTEGRATION  STEPS  IN  TIME  DOMAIN 
NXIMP:  NUMBER  OF  SPACE  STEPS  IN  IMPACT  RADIUS 
NTIMP:  NUMBER  OF  TIME  STEPS  IN  IMPACT  TIME 


2.  WITH  THE  ABOVE  DATA  FOLLOWING  PRIMARY  DATA  ARE  CALCULATED 
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RELEASE  2.0 


MAIN 


DATE  =  77139  21/11/39 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

C  3. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C  4. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C  5. 

c 

c 

c 

c 

c 

c 


B:  A  HALF  OF  THE  LAYER  THICKNESS 

K:  WAVE  NUMBER  FOR  FOURIER  TRANSFORM 
S:  LAPLACE  TRANSFORM  VARIABLE 

KO:  LIMIT  OF  INTEGRATION  FOR  INVERSE  FOURIER  TRANSFORM  FOR  X 
OMEGAO:  LIMIT  0<^  INTEGRATION  FOR  INVERSE  TRANSFORM  IN  TIME 
CO:  LAPLACE  TRANSFORM  PARAMETER 


DISPLACEMENT  AND  STRESS  FIELDS  ARE  CALCULATED  IN  TRANSFORMED  SPACE 
AND  BY  INVERSIONS  THESE  BECOME  DISPLACEMENTS  AND  STRESSES. 

THEY  ARE  GIVEN  IN  A  MATIX  FORM  AS  XX(I,J)  REPRESENTING  QUANTITY  AT 
ITH  TIME  STEP  AND  JTH  SPACE  STEP  IN  X 

U(I,J):  HRIZONTAL  DISPLACEMENT 
V(I,J):  VERTICAL  DISPLACEMENT 
SIGMAdtJJ:  NORMAL  STRESS 
TAU(I,J);  SHEAR  STRESS 
SIGMAKItJ):  TANGENTIAL  NORMAL  STRESS 


FOLLOWINGS  ARE  WORKING  MATRICES  FOR  THIS  PROGRAM 


DATAdfJ),  SUBd,J):  WORKING  MATRICES  FOR  SUBROUTINE  FOURT 
QO(I,J):  INTEGRAL  TRANSFORM  OF  IMPACT  FUCTION  GIVEN  IN  EQ(II-22) 

CAXdtJ),  CBXdtJ);  COS(ALPHA)  AND  COS(BETA)  IN  EQdI-16) 
XlBXd,J),  .YlAXd  ,J)  ,  .  .  :  Xl(BETA),  YKALPHA),..  IN  EQdI-18) 


FOLLOWING  SUBROUTINE  ARE  SUPPLIED  IN  THE  PRESENT  PROGLAM 

DPHASE:  CALCULATES  COS(BETA)  AND  CCS(ALPHA)  IN  EQ(  11-16)  WITH  GIV’ 
VALUES  OF  WAVE  NUMBER  K  AND  LAPLACE  TRANSFORM  VARIABLE  S 

DELL:  CALCULATES  XKEETA),  YKALPHA),.. 


IN  EQ  dI-19,20) 


oooooooooooooo 
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DET:  CALCULATE  0,D1,..  IN  £Q{II-23) 

FLIP:  THIS  PROGLAM  GENERATES  ONLY  A  HALF  OF  THE  DATA  (X>0)  DUE 
TO  SYMMETRY  AND  FLIP  THEM  TO  FIND  FULL  DATA  ALONG  X  AXIS 

MAP:  CONTROLLS  THE  PRINTOUT  FORMAT 

FACT:  TAKES  CARES  OF  THE  COORDINATE  SHIFT  AND  MULTIPLICATION  FACTORS 
IN  FOURIER-LAPLACE  DOUBLE  INTEGRAL  TRANSFORMS 

FOURT:  IS  FAST  FOURIER  TRANSFORM  ROUTINE  SUPPLIED  BY  IBM 


Q  a[c  ^  5j:  2}:  5;c  S';  3j:  ?!:  ^  ^ 

c 

c 


issl5:{c5 


>  ^  «lU 


:{c  5|;  sjc s{c  sj;  s^c  5»;  5{s  3|;  sjc  3}:  2|c  s^c  sj;  sjs  3{c  s^i  3{i Ji;  sjc  sjc  s 


IMPLICIT  REAL'T8(K) 

COMPLEX  DATA { 32 t 64) ,SUB(32t32) ,C AX( 32 » 32 ) » CBX ( 32 » 32 ) » QO ( 32 t32) 
COMPLEX  SIGMA{32,32) ,SIGMAl(32t32) ,TAU( 32»32) 

COMPLEX  UI32,32) ,V( 32,32 ) tVX( 32 t32) 

COMPLEX  Y1AX(32,32) , Y2 AX{ 32 , 32 ) , Y3AX ( 32 ,32 ) 

COMPLEX  X13X(32,32) , X2BX( 32 ,32 ) , X3BX ( 32 ,32 ) 

DIMENSION  NN(2),MM(32) 


C 

REALMS  C11,C12,C22,C66, CHAT, E66, V66,PI ,P2 
REAL=!=8  DCOS,DSIN,DBLE,DSQRT,DFLOAT,DLOG 

REALMS  DELTA ,D,RHO,FN, F,KX (32) , B , T AUO , A , DX , DT , UN  I TT, UN  I TX 
REALMS  Q,OMEGAO,BL,CO 
C 

C0MPLEX*16  CDEXP, COLOG, CDSQRTtCDCOS, CD  SIN 
C0MPLEX*16  BETA,ALPHA,CB,CA, SB,SA,C2NB,C2NA, S2NB, S2NA 
C0MPLEX*16  S,SI ,S2, 01,02,03, D4 
COMPLEX* 16  YIA, Y2A,Y3A,X1B,X2B,X3B 


C 

COMMON  YlA,Y2A,Y3A,XlB,X2B,X3e,Dl,D2,D3,D4,Cll,C12,C22,C66,CHAT 
EQUIVALENCE  (DELTA, D) 

EQUIVALENCE  ( DATA ( 1 , 1 ), SUB ( 1 , 1 ) ) 

C 

SI=(O.D  00,1.D  00) 

PI=3.1415926536D  00 
P2=PI*2.D  00 
INDEX=0 
INDEX=1 


C 

WRITE(6,4) 

4  FORMAT! ' 1* ///////////////20X ,' ***  WAVE  PROPAGATION  IN  COMPOSITE  P 

ILATE  ***• ) 
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5 

C 

C 

C 

C 

c 

c 

c 

c 

c 


WRITE(6,5) 

'=QRMAT(22X, 'GRAPHITE  F  I  BER  (  5  5^  )- EPOXY  MATRIX 


*  sis  5i!  #  sis  sic  sis  3l!  Jis  sj!  si!  jlt  sjt  *  * *  S|£  #  *  *  sjs  *  sic :!!  *  s|=  *  #  5tS  ❖  ^  sis  5ls  5r  Sr  *  '.S  ’I' - 


COMPOSITE* ) 

sls5!cs!ssjss|sslss!ssiss|cs(:#sls5ls 


INPUT  DATA  FOR  ELASTIC  PROPERTIES  OF  COMPOSITE  PLATE 

ALL  THE  DATA  ARE  SUPPLIED  IN  PSI  UNIT  BUT  NORMALIZED  BY  C66 

WHICH  IS  CONSTANT  REGARDLESS  THE  LAYUP  ANGLE 


sis  sis  sis  sis  :^s  sis 


-S  #  sis  sis  !{S  sic  sis  sic  sis 


READ  (5t101)  IANGLETClltCl2,C22»C66 
101  FORMAT! 110, 4D15. 7) 


CHAT=C11-C12*«2/C22 
IF  ( lANGLE.EQ. 100)  GO  TO  200 
WRITE (6,102)  I ANGLE,C11,C12,C22,C66 
102  FORMAT!//  20 X , • LA YUP  ANGLE  = ' , 1 3 , 3X , • DEGR EE  * 

$  /20X, 'Cd,  1)  =  ' ,D12.5,  '  PSI  •  ,10X, 'C!  1,2)  =  ' ,D12.5,  •  PSI' 

$  /20X, 'C(2,2)  =  '  ,D12.5, '  PS  I ' , lOX , • C ( 6 , 6 )  =  ' , D 12 . 5 , •  PSI'//) 

GO  TO  201 

200  CONTINUE 

WRITE{6,210)  Cll,Ci2,C22, C66 

210  F0RMAT(/20X,  '  PLATE  IS  ISOTROPIC  WITH  POISSCN"S  RATIO  1/4' 

$  /20X, 'C(l,l)  =  ' ,D12.5,  •  PSI • ,10X, 'C(l,2)  =  ' ,D12.5, •  PSI' 

$  /20X, *0(2,2)  =  ' ,D12.5, '  PS  I ' , lOX , ' C ( 6 , 6 )  =  • , D 12 . 5 , •  PSI'//) 

201  CONTINUE 
C 

E66='C66=!'6892.2D  00 

C11=C11/C66 

C12=C12/C66 

C22=C22/C66 

CHAT=CHAT/C66 

C66=1.D  00 


C 

Csicsississis^sis  WITH  CORRECTION  FACTOR 
C 

C66=P I**2/12.D  00 
C22=C22*C66 


C 

C 

C  - 

c 

C  INPUT  DATA  FOR  GEOMETRY  OF  COMPOSITE  PLATE 

C  ALL  THE  DATA  ARE  FIRST  SUPPLIED  IN  CGS  UNIT  BUT  CONVERTED  INTO  MKS  UNIT 

C 

C 

REA0(5,120)  NLAYER, DELTA, RHO 
120  FORMAT! 110,2020.10) 
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NLl=NLAYER+l 

FN=DFLOAT(NLAYER) 

8=DELTA/FN 

WRITE (6, 121 )  DELTAtRHO,NLAYER,e 

121  F0RMAT(20X, 'TOTAL  THICKNESS  OF  COMPOSITE  PLATE  ;  DELT A= • , F 10 . 5 , 

$  '  CM'/ 

$  20X, 'DENSITY  OF  COMPOSITE  ;  RHO=  '  ,  F 10. 5,  '  GR/CM*=!=3'/ 

$  20X,  'PLATE  IS  MADE  OF  '»I3,3X,  'IDENTICAL  LAYERS'/ 

$  20X, 'LAYER  THICKNESS  ;  2B=',F10.5,'  CM'//) 

C 

DELTA=DELTA/100.D  00 
RH0  =  RH0>i'1000.D  00 
B=B/200.D  00 
C 

c 

C  - - - 

c 

C  INPUT  DATA  FOR  IMPACT 
C 

c 

READ(5,60)  NA,TAUO 

60  FORMAT!! 10,D20.10) 

READ(5,111)  NXtNT,NXIMP,NTIMP 

111  F0RMAT(4I10) 

C 

WRITE{6,112)  NX,NXIMP»NT,NTIMP 

112  F0RMAT(20X,' TOTAL  SPACE  STEPS;  NX  = ' , 1 3 , 5X , ' W ITH' , 13 , 2X , ' STEPS  FOR 
$  CONTACT  RADIUS'/ 

$20X, 'TOTAL  TIME  STEPS  ;  NT= • , 13 , 5X , ' W ITH ' , 1 3 ,2X , ' ST EP S  FORCONTAC 
$T  TIME'//) 

C 

A=DFLOAT(NA)*DELTA 

DX=A/DFLOAT(NXIMP) 

DT=TAUO/DFLOAT(NTIMP) 

C 

WRITE(6,61)  A,DX,TAUOtDT 

61  F0RMAT(20X, 'CONTACT  RADIUS  ;  A=',D12.5,'  M'/ 

$20X, 'SPACE  STEP  ;  DX=',D12.5,'  M'/ 

$20X, ' CONTACT  TIME  ;  T AUO= ' t D 12 .5 , •  SECOND'/ 

$20X,'TIME  STEP  ;  DT=',D12.5,'  SECOND'//) 

C 

C 

C  - 

C 

C  NORMALIZE  ALL  THE  INPUT  DATA 

C 

C 


V66=DSQRT( E66/RH0) 


ooooooooooo 
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UNITT=A/V66 

UNITX=DELTA 

F=D*D*RH0/ (E66*UNITT**2)/2.D  00/FN 

NX2=NX/2 

NN(1)=NT 

NN(2)=NX 

DX=DX/UNITX 

0T=DT/UNITT 

K0=PI/DX 

OMEGAO=PI/DT 

BL=A/D 

C0=DL0G{1.D  06*2. D  00*DT)/(3.D  00*DT*DFL0AT( NT) ) 


-r*  ■*r 't'  ■T'  ■n*  'r*  "f*  ^  ^  -r  'r*  'T*  'r  "r  'i'  A* 


^  wu  ^ 
^  ^ 


j)j  jAj  jlj  jjj  jl* 


V* 


r  'I*  *Y*  -r* 


^  'f-  .T(»  -p  ^  ^ 


CALCULATE  THE  IMPACT  INPUT  FUNCTION  QO(I,J)  IN  EQI 11-22) 

CALCULATE  COS(6ETA)  AND  COS(ALPHA)  IN  EQI 11-16)  BY  SUBROUTINE  DPHASE 
CALCULATE  XKBETA),  YKALPHA)  ,..  BY  SUBROUTINE  DELL 


DO  30  J=1,NX2 

K=2.D  00*K0*(DFL0AT(J )-.5)/DFL0AT(NX)-K0 

K2=K**2 

KX( J)=K 

Q=PI**2/DSQRT( P2)*DSIN(K*BL)/K/( ( K*BL ) **2-P 1**2 ) 

C 

DO  30  I=ltNT 

S=C0+SI*0MEGA0*(1.D  00- ( DFLOAT ( I )-.5D  00)*2.D  OO/DFLO AT{ NT ) ) 
S2=S**2*F 

Q0(ItJ)=Q/2./S*(l.D  OO-CDEXP I-S*TAUO/UNITT) ) * ( P2*UNI TT ) **2 
$  /( {S*TAU0)**2+(P2*UNITT)**2) 

C 

CALL  DPHASE  ( K , S2 , CB , C A , NLAY ER ) 

CBX( I , J)=C8 
CAXI I , J)=CA 
C 

CALL  DELL(K,S2,CB,CA,SI,NLAYEP ) 

Y1AX( I , J)=Y1A 
Y2AX( I tJ)=Y2A 
Y3AX( I , J)=Y3A 
XIBXI I ,J)=XLB 
X2BX( I , J)=X2B 
X3BX( I , J)=X3B 
C 
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30  CONTINUE 
C 

c 

c 

£  3!j  ^  Jit  5|t  j;;  Jit  :is  ^  5|S  Jl;::!!*  5^5  Jjs  51;  >|£  s|;  5}:  sj:  ^  :is5!s  Jjc  :is  :j!  j!:*  5|!  5|s  =|:  5!:  51:  ^  =!: 

c 

c 

C  REPRODUCE  THE  IMPACT  FUNCTION  TO  CHECK  INPUT 
C 

DO  300  1=1, NT 
DO  300  J=1,NX2 

300  SUB( I,J)=QO( I,J) 

CALL  FLIP(DATA,NX,NX2,NT,+1) 

CALL  F0URT(DATA,NN,2,-1,1,0) 

CALL  FACT(DATA,NX,NT,CO,OMEGAO,KO,PI,SI) 

WRITE(6,301) 

301  FORMAT  (•  l'///////////////20X,  REPRODUCTION  OF  IMPACT  FUNCTIO 

IN  #««») 

CALL  MAP(DATA,NX,NT,NX2 , MM, INDEX) 

C 

C 

c 

c** 
c 
c 
c 
c 
c 
c 
c 
c 
c 


5;s?*;2»::{C5;s3;;5{c:{s  3{s:}c  2^  3!^:^C3;^3^£5!^  :?s:js5jt35ss;r  5*c:*s5;c5js 


THIS  IS  THE  MAIN  PART  OF  THE  PROGRAM. 

CALCULATE  D(  BETA)  ,OBAR(  ALPHA)  .  IN  EQni-19,20)  BY  SUBROUTINE  DET 
NEXT  CALCULATE  1,V,..  IN  EQ(II-18)  IN  TRANSFORMED  SPACE 
AND  FLIP  TO  FIND  FULL  DATA  AND  INVERT  THEM  BY  MEANS  OF  FOURT. 
REPEAT  THIS  PROCESS  FROM  N=0  TO  NLAYER 

DO  11  N=1,NL1 

NY=N-1 

NYY=NY-1 


C 

C 

C  - 

C 

C  GENERATION  OF  DATA  FOR  DISPLACEMENTS  AND  STRESS  IN  TRANSFORMED  SPACE 
C 

DO  100  J=1,NX2 
DO  100  1=1, NT 
CB=CBX(  I  ,J) 

CA=CAX(  I  ,J) 

X1B=X1BX(I,J) 

X28  =  X2eXI I  ,J  ) 

X3B=X3BX(  I  ,J  ) 
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Y1A  =  Y1AX( I  ,  J) 

Y2A=Y2AX(  I  ,  J) 

Y3A=Y3AX( I , J) 

C 

SB=C0SQPT( l.D  00-CB**2) 

SA=CDSGRT( l.D  00-CA**2) 

BETA=CB+SI*SB 
ALPHA=CA+SI*SA 
BETA=CDLOG{BETA)/SI 
ALPHA =CDLOG( ALPHA) /SI 
C 

CALL  OET(ALPHA,BETA»SI,FN) 

C2N8  =  CDC0S(2.D  00*8ETA>!=DFLOAT(NY)  ) 

S2N8  =  CDSIN(2.D  00*BETA=!=DFL0AT(  NY)  ) 

C2NA  =  CDCOS(2.D  00=^ALPHA«DFLCAT  ( NY  )  ) 

S2NA  =  CDSIN(2.D  00=!'ALPHA*DFL0AT  (NY  )  ) 

C 

U(  I,  J  )  =  (  XlB*(Dl*C2NB  +  SI*D2=!=S2NB)+YlA*(04=f=C2NA+SI=!=03*S2NA)  )*Q0(  I  ,J) 
V{ I , J) =(X2B*(D2*C2NB+SI*D1*S2NB)+Y2A#(D3*C2NA+SI«D4*S2NA) )«Q0( I , J) 
TAU{  I  tJ)=(  X38*{D2*C2NB  +  SI*D1*S2NB)  +  (03*C2NA+SI*D4«S2NA)  )=f=Q0(  I,  J) 
SIGMA(  I  »J)=(  (D1*C2N8+SI*D2*S2NB)+Y3A*(D4*C2NA+SI*D3’«*S2NA)  )*Q0(  I ,  J  ) 
IGO  CONTINUE 
C 
C 

C  - 

c 

C  INVERSION  AND  PRINTOUT  OF  HORIZONTAL  DISPLACEMENT  UN(I,J) 

C 

DO  10  1=1, NT 
DO  10  J=1.NX2 
10  SUB( I , J)=U( I  ,J  ) 

CALL  FLIP(0ATA,NX,NX2,NT,-1) 

CALL  FOURT(DATA,NN,2,-1,1,0) 

CALL  FACT(DATA,NX,NT,CO,OMEGAO,KO,PI  ,SI  ) 

WRITE(6,981)  NY 

981  FORMAT!* l'///////////////20X,'U' ,13) 

CALL  MAP(DATA,NX,NT,NX2,MM, INDEX) 

C 

C 

C  - 

c 

C  INVERSION  AND  PRINTOUT  OF  VERTICAL  DISPLACEMENT  VN(I,J) 

C 

DO  20  1=1, NT 
DO  20  J=1,NX2 
20  SUB( I , J)=V( I ,  J  ) 

CALL  FLIP(DATA,NX ,NX2,NT,+1) 

CALL  FOURT(DATA,NN,2,-1,1,0) 


-  85  - 


1  RELEASE  2.0  MAIN  DATE  =  77139  21/11/39 

CALL  FACTIDATA,NX,NT,C0,0MEGA0,K0,PI»SI  ) 
k^RITE{6,982)  NY 

982  FORMAT! ' 1* /////// //////// 20X V  t 13) 

CALL  MAP(DATAtNX,NT,NX2t MM, INDEX) 

C 

C 

Q  - 

C 

C  INVERSION  AND  PRINTOUT  OF  SHEAR  STRESS  TAU(I,J) 

C 

DO  35  1=1, NT 
DO  35  J=1,NX2 
35  SUB( I , J)=TAU( I , J) 

CALL  FLIP(DATA,NX,NX2,NT,-1) 

CALL  FOURT(DATA,NN,2,-1,1,0) 

CALL  FACT{DATA,NX,NT,CO,OMEGAO,KO,PI,SI ) 

WRITE(6,983)  NY 

983  FORMAT (• 1* ///////// //////20X ,' TAU' ,13) 

CALL  MAP(DATA,NX,NT,NX2,MM, INDEX) 

C 

C 

c  - - - 

c 

C  INVERSION  AND  PRINTOUT  OF  NORMAL  STRESS  SIGMA(I,J) 

C 

DO  40  1=1, NT 
DO  40  J=1,NX2 
40  SUB{ I , J)=SIGMA( I,J) 

CALL  FLIP(DATA,NX,NX2,NT,+1) 

CALL  FCURT(DATA,NN,2,-1,1,0) 

CALL  FACT! DATA, NX,NT,CO,OMEGAO,KO,PI ,S  I ) 

WRITE{6,984)  NY 

984  FORMAT! ' 1' ///////////////20X , ' SIGMA' ,13) 

CALL  MAP  ID AT A, NX, NT  ,NX 2, MM, INDEX) 

C 

c 

Q  - 

c 

C  INVERSION  AND  PRINTOUT  OF  TANGENTIAL  NORMAL  STRESS  SIGMA1!I,J) 

C 

IP  INY.EQ.O)  GO  TO  160 
DO  50  1=1, NT 
DO  50  J=1,NX2 

50  SUB!  I  ,  J)=SIGMAU  I  ,  J  )  +  FN*C  1 2*  V(  I,  J) 

CALL  FLIP!DATA,NX,NX2,NT,+1) 

CALL  F0URT!DATA,NN,2,-1,1,0) 

CALL  FACT!DATA,NX,NT,CO,OMEGAO,KO,PI,SI ) 

WRITE!6,985)  NYY 


RELEASE  2.0 


MAIN 


DATE  =  77139 


21/11/39 


985  FORMAT! • 1* ///////////////20X,' SIGMAl' , 13) 

CALL  MAP (DATA, NX, NT  ,NX2 , MM , I NDEX ) 

C 

IF  ( NY.EQ.NLAYER)  GO  TO  70 
DO  51  1=1, NT 
DO  51  J=1,NX2 

51  SIGMAl (I ,J)=-SI*KX( J)*C11*U( I,J)-FN*C12#VX( I, J) 

GO  TO  80 
C 

160  CONTINUE 

DO  161  1=1, NT 
DO  161  J=1,NT 

161  SIGMAK  I  ,J)=-SI*KX(  J)*C11*U(  I,J)-FN#C12«V(I,J) 

GO  TO  80 

C 

70  CONTINUE 

DO  71  1=1, NT 
DO  71  J=1,NX2 

71  SUB(  I  ,  J)  =-SI*KX(J  )=!'C11*U(  I  ,  J)  +  FN*C12*(  V!  I  ,  J)-VX(  I  ,J)  ) 
CALL  FLIP(DATA,NX,NX2,NT,+1) 

CALL  F0URT(DATA,NN,2,-1,1,0) 

CALL  FACT(DATA,NX,NT,CO,OMEGAO,KO,PI,SI) 

WRITE(6,985)  NY 

CALL  MAP(DATA,NX,NT,NX2,MM, INDEX) 

GO  TO  90 
C 

80  CONTINUE 

DO  81  1=1, NT 
DO  81  J=1,NX2 

81  VX( I , J)=VI I , J) 

C 

90  CONTINUE 

C 

c 

11  CONTINUE 

C 


-T*  T*  -r  ■n*  ^  •O'  'r»  •T' 


ajc  aj;  3}:  51:  ^  *  3^  5^  3{:  ^ 


'  'j'  ^  'r*  -r  'p  O'  'r  'T*  o'  '■r  -r  o- 


C 

c 

c 

c 


^  tJL. 
^  ^ 


STOP 

END 
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sf:  *  5}:  *  *  3!:  sis  *  5|5  *  ’f' 


sis  sftslssjssls^jsslsslei 


s  S|5  sjs  *  sjs  s|s  sjs  sis  s|c  $  si:  #  sic  sf:  sis  :{s  s|s  sis  sis  ^ 


THIS  SUBROUTINE  CALCULATES  THE  PHASE  SHIFT  BETA  AND  ALPHA  FROM 
EQ( 11-15, 16)  OF  THE  PRESENT  REPORT  WITH  GIVEN  VALUES  OF 
WAVE  NUMBER  K  AND  LAPLACE  TRANSFORM  VARIABLE  S 


CA=COS( ALPHA) 
CB=COS(BETA) 


w  -UU 


>  'JU  .A.  JU  .1#  ^  ^  .A.  « 

.  ^  ^  ^  iff*,  ^  ^ 


^  -v  'r  ‘ 


^  'T'  'Y'  n'  nr  »r‘ 


JU  ^ 

nr  "r  ^ 


nr  ^  nr  nr  "r 


n'  ■•'t'  n"*  n*  n*  n*  •'r  n'*  n*  ^  n^  n*  nr  nr  n'  n*'  nr  •'i'  •'r  nr  nr 


SUBROUTINE  DPHASE(K,S  , CB , C A,NLAY EP ) 

IMPLICIT  C0MPLEX*16{A,X,Y) 

C0MPLEX*16  R00TP,R00TM,S,CDSQRT,DCMPLX  ,DCB,DCA 
C0MPLEX*16  D1,D2,D3,D4 
C0MPLEX*16  CB,CA 

REALMS  C11,C12,C22,C66,CHAT,N,K2,DFL0AT ,DBLE,K 
C 

COMMON  Y1A,Y2A,Y3A,X1B,X2B,X3B,D1,D2,D3,D4,C11,C12,C22,C66,CHAT 
ROOTP(AA,AB,AC)  =  {-AB+CDSQRT(  AE=!=*2-4.D  00’!=AA=i=AC)  )/ (  2.D  00*AA) 
ROOTM(  AA,  AB,  AC)  =(-AB-CDSQRT(  AB**2-4.D  00*AA*AC )  )  /  (  2. D  OOT-AA) 

C 

N=DFL0AT(NLAYER)*2.D  00 
K2=K**2 
C 
C 

A1=(K2*C11/N+S)*(K2*C66/N+S) 

A2=K2*CHAT/( 3.D  00-N)+N*C66+S/3.D  00-C 12 *C66*K2/ ( 3 . D  00*N-C22) 
A2  =  A2*(-C12*K2/(3.0  00=!'N)+N*C22  +  S/3.D  00) 

A3=N*C22+S/3  .D  00-K  2-C12-  (  C66s|sK2/N+S  )  /  (  9 .  D  00*N**2=!'C2  2  ) +C66*K2/ 
$  (3.D  00*N) 

A3  =  A3sis(Cll=4=K2/N  +  S)-K2*(  C12  +  C66)**2 

A3=A3  +  (C66*K2/N+S)’i={CHAT=isK2/(3  .D  00*N)  +N’i^C66  +  S/3 . +C  12 **2*K2/ 

$  (3.D  OO^CZZ)) 

C 

AA=A1+A2-A3 
A8=A3-2.D  00*A2 
AC=A2 

DCB=ROOTPI AA , AB,AC) 

DCA=ROOTM( AA,AB,AC) 

C 

CB=CDSQRT(DCB) 

CA=CDSQRT(DCA) 

RETURN 

END 


OOOOOOOOOOOOOO 
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iJU  .JU  <JU  UU  JL.  'JU  iJU  ^  UU  Ui.  -JU 


.A#  4jl^  vl^  .A.  V.-  .A.  .A^ 


^1^  %!*.  ^A. 

-r  *>'•  -n'  -t'  'f  '1'*  "r 


.A^  JU  .A^ 


«L.gL><4>u«gL>JU<a«a«<jU'./ua.^>A.<a^sk 
n'  '!'•  *T*  n*  -nr  n-  nr  n'  'rnr  ■»('  n' 


THIS  SUBROUTINE  CALCULATES  DELTA ( BETA) , DELTABAR { ALPHA ), DELTAl ( BETA ),. . 
IN  EQ{II-19t20)  and  XI { BETA ) ,Y 1 { ALPHA )  IN  EQ ( 1 1-18 , 19 , 2 0) 

X1B=X1(BETA) ,  Y1A=Y1( ALPHA) 


Uf  %l«»  OU  WU  WU  ^  JU  WU  ^1#  V-*  %W  VU  WU  WU  WU  WU  «l.i*  %iU  WU  4JU  <ijl^  JU  ^jU  WU  WU  WU  WU  «Ml«r  WU  JU  0«r  >IU  wu  <ou  ^ 


'  5r  ^ 


SUBROUTINE  DELL  ( K , S , CB ,C A , S I , NLAY ER ) 

IMPLICIT  REAL*8(K) , COMPLEX* 16 ( S , X , D» Y ) 

C0MPLEX*16  SB,CB,CDSQRT,SA,CA 
REAL*8  NtDFLOAT 
REAL*8  C11tC12tC22,C66,CHAT 
C 

COMMON  Y1A,Y2A,Y3A,X1BtX2B,X3E,D1,D2,D3,D4,C11,C12» C22,C66,CHAT 

c 

K2=K**2 

N=OFLOAT{NLAYER) 

C 

C 

S11=S+C11*K2/2.D  00/N 
S66=S+C66*K2/2.D  00/N 
SA=CDSQRT(1.D  00-CA**2) 

SB=CDSQRT( l.D  00-CB**2) 

C 

C 

DELTAB=C6**3*S11*S66+SB**2*CB*(-C66*(C66+C12)*K2 
$  +S66*(S/3.D  00+CHAT*K2/6.D  OO/N+2.0  00*N*C66)) 

DELie=SI*K*SB**2*CB*(-C66-C12+C12*S66/6.D  00/N/C22 ) 
DEL2B=SI*SB**3*( C12*C66*K2/6.D  00/N/C22-S/3 . D  00 
$  -CHAT*K2/6.D  00/N-2.D  00*N*C66 )-S I *CB**2*S B*S 1 1 
0EL3B=C8**2*SB*(C12*K*S11*S66/6.D  00 /N/C22-C 66*K*S 1 1 ) 

$  +SB**3*(C12*K*(S/3.D  00+CHAT*K2/6.D  00/N+2.D  00*N*C66) 

$  -C12**2*C66*K*K2/6.D  00/N/C22) 

C 

0ELTAA  =  SI*CA**2*SA*{ ( C12+C66)*C66*K2-S66*( S/3  .D  00+CHAT*K2 
$  /6.D  OO/N+2.0  00*N*C66+C12**2*K2/6.D  00/N/C22)) 

$  +SI *SA**3*(S/3.D  00+2. D  00*N*C22 ) * ( C66 *C 1 2*K2/6 . D  00/N/C22 
$  -S/3.D  00-CHAT*K2/6.D  00/N-2.D  00*N*C66) 

DEL1A=SA**2*CA*(-K2*C12*S66/6.D  00/6. D  00/N**2/C22+C6 6*K2 /6. D  OO/i; 
$  +S/3.D  00+2. D  00*N*C22)+CA**3*S66 

DEL2A=CA**2*SA*K*(C12+C66) 
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$  +SA**3*{K/6.D  OO/N* ( S/3.D  00+CHAT*K2/6 . D  00/M+2.D  00*N*C66) 

$  -C66*C12*K**3/36.D  00/C22/N**2} 

DEL3A=-S  I*CA«*3*C12*K*S66+S  I*SA=«'*2’!=CA’!'(  C66*K=!=  ( S/3  .D  0  0+2.D  00*N 
$  =!=C22+C66*K2/6.D  00/N)-K/6.D  00/N=!=S66T-{  S/3  .D  OO+CHA  T*K2/6.  D  00 

$  /N+2.D  00*N*C66) ) 

C 

C 

X1B=-DEL1B/DELTAB 

X2B=-DEL2B/DELTAB 

X3B=-DEL3B/DELTAB 

Y1A=-DEL1A/DELTAA 

Y2A=-DEL2A/DELTAA 

Y3A=-DEL3A/DELTAA 

C 

RETURN 

END 


o  o  o  O  O  o 
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*1,.  ^  .jU  ■Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  mJLr  Jr  Jr  Jr  -Jr  Jr  Uir  Jr  Jr  Jr  Jr  ■Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  Jr  ■Jr  Jr  jl.  ^  Jr  JJ  ^  r 

rx^  r^^rf,  r^^rf.  .y'  ^  'T'  n'  ^  *1'  ^  .y.  ^  'I'  'I'  '1^  'T»  ^T*  'T*  ^  .1,*  .y.  ^  ^  ^  ^  ^  ^  rf.  rf.  r^-,  ^,x  r^-  rf^  r 


•  'r-  -n'  “T*  -r  -r 


xU  -JU  xJ 

rf*  rf*  r 


THIS  SUBROUTINE  CALCUALTES  D»01,..  IN  EQ( 11-23)  OF  THE  PRESENT  REPORT 

£  sit  Jis  :(s  >|s  5|S  :!s  :it  :j!  *  3|J  J|S  ;!j  $ :(!  :!c  #  «  s|5  4:  *  5|«  *  *  !fs  sf:  i(5  #  a!*  XS  =r  *  #  Sr  *  #  5r  *  #  *  Sc  *  ^  s|s  *  s|=  #  Sr  *  5|s  # ’!' >!«  Sr  #  *  Sr  *  *  *  Sr  »|s  #  *  =(S * 

C 

C 

SUBROUTINE  DET ( ALPH A , B ETA , S I » FN ) 

IMPLICIT  C0MPLEX#16(DtX,Y) 

C0MPLEX*16  ALPHA, BETA 

COMPLEX# 16  C2Ne,C2NA,S2NA,S2NB,CDSQRT,SI 
REAL#8  FN 

REAL#8  C11,C12,C22,C66,CHAT 
C 

COMMON  Y1A,Y2A,Y3A,X1B,X2B,X3B,D1,02,D3,D4,C11,C12,C22,C66,CHAT 
C 

C2NA=CDCOS(2.D  00#ALPHA#FN) 

S2NA=CDSIN(2.D  00*ALPHA#FN) 

C2NB=CDC0S(2.D  00#BETA#FN) 

S2NB=CDSIN(2.D  00#BETA#FN) 

C 

X=Y3A*X36#(1.D  00-C2NAsisC2NB ) 

Y=X3B#Y3A#S2NB#C2NA-S2NA#C2NB 

C 

D=-2.D  00#X+(1.D  00+X3B##2#Y3A##2)#S2NA#S2NB 
D1=-(X-S2NA*S2NB) 

D2=-SI#Y 

D3=SI#Y*X3B 

D4=X3B#( X3B#Y3A#S2N6#S2NA+C2NA#C2NB-1.D  00) 

C 

D1=D1/D 

D2-D2/D 

03=03/0 

04=04/0 

C 

RETURN 

END 


oooooooooo 
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C 

C 


ALL  THE  DATA  IN  THE  MAIN  PROGRAM  ARE  GENERATED  FOR  ONLY  HALF  OF  THE 
PLATE  WHEN  X>0.  DUE  TO  SYMMETRY  OF  THE  PROBLEM  WE  CAN  GENERATE 
THE  FULL  DATA  BY  FLIPPING  THE  HALF  OF  THE  DATA. 


:5c sjc 5|s 5{i  sits;: :^rs 


1 3tc  sj: :^s 3jc  35: 3?:  :5s  3{c  sit  35:  3}:  sj:  :ic  :Cc 5*;  :^t  sjt  sjt  s5t  sfc  ^ 


^  JU  *iU  <.1^ 


,1,  ^  ^■JL. 


<Jj.  yir  « 

•'t*  -v  •nr  'c  ' 


SUBROUTINE  F L I P ( DAT A , NX , NX2 t NT , I NDEX ) 
COMPLEX  DATA(NT,NX) 

DO  10  J=1,NX2 
JJ=NX+1-J 
DO  10  I=L»NT 

DATA! I,JJ)=FLOAT( I NDEX ) *DAT A ( I  ,  J ) 

10  CONTINUE 
RETURN 
END 
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C 

C 

Q  ««  sjs  !|;  *  sis  *****:<«  5(5  *  ^  :!s  J|s  sis  Ji!  jj!  :j!  **  :i5  #  St:  :ic  #  :j:  sjc  #  :jt  5|s  sis  sU  «  5|s  :!£  Jjc  sjc  :(s  *  sis  *  Jls*#**  sj!  * 

C 

C 

C  THIS  SUBROUTINE  TAKES  CARE  OF  THE  COCRIOINATE  SHIFT  IN 
C  LAPLACE  AND  FOURIER  TRANSFORM  IN  THE  PROCESS  OF  APPLYING 

C  FAST  FOURIER  TRANSFORM  ALGORITHM 

C 
C 

:<t  sis  :(s  sfs  sis  3|s  s(s  s|s  #  sis  #  4s  sjs  sis  s|s  sis «  s|s  sjs  *  sis sjs  sjs  sjs  sjs  s|s  sis  *  sis  sis  :t:  *  :js  sis  s^  :{s  sjs  St:  5|s  ;(s  sjs :{:  *  3|s  !js  s|s  3(s  sjc  sis  s|s  s|5  *  3(s  sjc  s|t  s|s  :js  s|t  sjs « sis  sjs  :{s  sis  s|s  :ts  sis  sjs  sis  sis 


C 

C 

SUBROUTINE  F ACT ( DAT A , NX ,NT , CO , WO , KO t P I » S I ) 
COMPLEX  DATA(NT,NX) 

C0MPLEX*16  CDEXPtSI 
REAL*8  DEXP,DSQRT,DFLOAT 
REAL*8  COtWOsPI ,KO,FT,FX 
FX=DFLOAT{NX) 

FT=DFLOAT(NT) 

NX2=NX/2 


C 

DO  10  L=1,NX2 
DO  10  M=1,NT 

DATA!  M,L)  =  DATA(M,L)*4.D  00*K0=!=W0/ {  DSQRT  (  2  .  D  00*PI  )  **3=!=FT=f=F  X) 
$  ^DEXPICO^PI^DFLOATIM-D/WOI^CDEXPISI^PK'!  l.D  00-1. D  00/FX) 

$  *0FL0ATIL-1)  )*CDEXP(SI’!=PI*(  l.D  00-1. D  00/FT)  *DFLOAT(  M-1)  ) 

10  CONTINUE 
RETURN 
END 
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C  sis  Jls  5|:  Jls Jl! « sj!  s|!  **  s|5  =1=  ifj  «***«  =!=«  sfs  >1=  s|s  *  «***«**’!==!«*****  ^  =1=  * 

C 

C 

C  THIS  SUBROUTINE  CONTROLS  THE  FORMAT  OF  THE  PRINTOUT  OF  THE  FINAL  RESULTS 
C 

C  IF  INDEX=0:  ALL  THE  NUMERICAL  VALUES  ARE  PRINTED 

C  =l:  THE  MAXIMUM  AND  NORMALIZED  VALUES  ARE  PRINTED 

C 
C 

{;;s|s*«s!s5!s*s|ss!s5|s$s!::js*:{ssls#s|ss«:s!s#s!ssis3!s#slss!s3lsslc3!s*:{5«8!s*slss|s«sJs««5issls«s!ss|s*s!s:Ss|ss!s3!s3|cslc5!ss!ss|ss|s5|s3iss|s5!ssjss|ss!ss!ss!s«:^ssiss|t:{s 


c 

c 

SUBROUTINE  MAP (DATA , NX, NT tNX2, MM, INDEX) 

COMPLEX  DATA(NT,NX)  ,S 
DIMENSION  MM(NX2) 

IF  ( INDEX. EQ.l)  GO  TO  200 
DO  44  IQ=1,NT 

44  WRITE(6,15)  I Q , ( D AT A ( I Q , I ) , I =1 ,NX2 ) 

15  FORMAT! I5,2E14.5,2X,2E14.5, 2X,2E14.5,2X,2E14.5/ 

$  3(5X,2E14. 5,2X,2E14.5,2X,2E14.5,2X,2E14.5/)/ 

$  4(5X,2E14.5,2X,2E14.5,2X,2E14.5,2X,2E14.5/)/) 

200  CONTINUE 

Csfssjssis  FIND  the  MAXIMUM  VALUE 
RS=  l.E-3 
NT5=NT-5 
NX5=NX2-5 
DO  114  1=1, NT5 
DO  114  J=  1,NX5 
S=  DATA! I, J) 

TP=  REAL!S)/RS 

IF  ! ABSI TP) .LT.l.  )  GO  TO  114 
RS=  REAL(S) 

114  CONTINUE 

WRITE  !6,516)  RS 

516  F0RMAT!20X,  •-**  MAXIMUM  VALUE  =',E12.5,‘ 

DO  119  1=1, NT5 
DO  113  J=1,NX5 
S=  DATA! I , J) 

113  MM( J)=  REAL! S)/RS*100 

WRITE!6,515)  ( MM ( K I M) , K IM= 1 , 27 ) 

119  CONTINUE 

515  FORMAT! 10Xt27I3) 

RETURN 

END 
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1 


2 

C 

C 

C 


C 

C 

C 

5 


10 

11 

12 


20 

30 


31 

32 


40 

50 

51 

60 

C 

C 

C 


subroutine  F0URT(DATA,NN,ND IM, is  I GN, I  form 

,WOPK) 

FFTTOOO 

DIMENSION  DATAI 1) ,NN( 1) , I  FAC T ( 32 ) , WORK ( 1 ) 

FFTT077. 

TW0PI=6. 283185307 

FFTT078 

IF(N0IM-1)920, 1,1 

FFTT079C 

NT0T=2 

FFTTOSOl 

DO  2  IDIM=1,NDIM 

FFTT081 . 

IF(NN( IDIM) ) 920,920,2 

FFTT082C 

NTOT=NTOT*NN( IDIM) 

FFTT083C 

FFTT084. 

MAIN  LOOP  FOR  EACH  DIMENSION 

FFTT085C 

FFTT086C 

NP1  =  2 

FFTT087i. 

DO  910  IDIM=1,N0IM 

FFTT088C 

N=NN{ IDIM) 

FFTT089" 

NP2=NP1#N 

FFTT090C 

IF(N-1)920,900,5 

FFTT091C 

FFTT092 

FACTOR  N 

FFTT093C 

FFTT094( 

M=N 

FFTT095  . 

NTW0=NP1 

FFTT096C 

IF  =  1 

FFTT097C 

IDIV=2 

FFTT098C 

IQUOT=M/IDIV 

FFTT099C 

IREM=M-IDIV4IQU0T 

FFTTIOOC 

IF( IQUGT-IDIV)50, 11,11 

FFTTIOLC 

IF( IREM)20,12,20 

FFTT102C 

NTWO=NTWO+NTWO 

FFTT103': 

M=IQUOT 

FFTT104C 

GO  TO  10 

FFTT105C 

I0IV=3 

FFTT106' 

IQUOT=M/IDIV 

FFTT107C 

IREM=M-IDIV*IQUOT 

FFTTlOet 

IF( IQUCT-I0IV)60,31,31 

FFTT109': 

IREM)4C,32,40 

FFTTllO 

IFACT(  IF)  =  IDIV 

FFTTlilc 

IF=IF+1 

FFTT112C 

M=IQUOT 

FFTT113-, 

GO  TO  30 

FFTT114C 

IDIV=IDI V+2 

FFTT115( 

GO  TO  30 

FFTT116C 

IF{ IREM)60,51,60 

FFTT117C 

NTWO=NTWO+NTWO 

FFTT118: 

GO  TO  70 

FFTT119C 

IFACT( IFJ=M 

FFTT120C 

FFTT121C 

SEPARATE  FOUR  CASES — 

FFTT122C 

1.  COMPLEX  TRANSFORM  OR  REAL  TRANSFORM 

FOR  THE  4TH,  5TH,ETC. 

FFTT123C 
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C 
C 
C 
C 

c 

70  N0N2  =  NP1’!=(  NP2/IMTW0) 

ICASE=1 

IF(IDIM-4)71,90,90 

71  IF( IF0RM)72,72,90 

72  ICASE=2 
IF(IDIM-1)73,73,90 

73  ICASE=3 

IF( NT WO-NP 1)90,90, 74 

74  ICASE=4 
NTW0=NTW0/2 
N=N/2 
NP2=NP2/2 
NT0T=NT0T/2 
1  =  3 

DO  80  J=2,NT0T 
DATA( J)=DATA( I ) 

80  1=1+2 

90  I1RNG=NP1 

IF(  ICASE-2) 100,95, 100 
95  I1RNG=NP0*( l+NPREV/2) 

C 

C  SHUFFLE  ON  THE  FACTORS  OF  TWO  IN  N 

C  CAN  BE  DONE  BY  SIMPLE  INTERCHANGE, 

C 

100  IF(NTWO-NP1)600,600,110 

110  NP2HF=MP2/2 

J=1 

DO  150  I2=1,NP2,N0N2 
IF( J-I2) 120, 130, 130 
120  IlMAX=I2+NON2-2 

DO  125  11=12, UMAX, 2 
DO  125  13=11 ,NT0T,NP2 
J3=J+I3-I2 
TEMPR=DATA( 13) 

TEMPI=DATA{  13+1  ) 


FFTT124C 
FFTT125C 
FFTT126C 
FFTT127C 
FFTT128C 
FFTT129C 
FFTT1300 
FFTT13'1C 
FFTT132C 
FFTT133G 
FFTT134C 
FFTT135C 
FFTT136C 
FFTT137C 
FFTT138C 
FFTT139C 
FFTT140C 
FFTT141C 
FFTT142C 
FFTT143C 
FFTT1440 
FFTT145C 
FFTT146C 
FFTT147<: 
FFTT148C 
FFTT1490 
FFTT150C 
FFTT151C 
FFTT152C 
FFTT153C 
FFTT154C 
FFTT155: 
FFTT156C 
FFTT157C 

AS  THE  SHUFFLING  FFTT158C 
NO  WORKING  ARRAY  IS  NEEDED  FFTT159C 

FFTT160: 
FFTT161C 
FFTT162C 
FFTT163C 
FFTT164C 
FFTT16.5C 
FFTT166C 
FFTT167C 
FFTT168C 
FFTT:.  69C 
FFTT170C 
FFTT171C 


C 

C 

C 

C 

C 

C 

c 

c 


THE  2ND  OR  3RD  DIMENSION. 
DATA,  SUPPLYING  THE  OTHER 


METHOD- 
HALF  BY  CON- 


DIMENSION, 
EACH  STAGE, 


N  ODD.  METHOD- 
SUPPLYING  THE  OTHER 


4. 


DIMENSIONS. 

REAL  TRANSFORM  FOR 
TRANSFORM  HALF  THE 
JUGATE  SYMMETRY. 

REAL  TRANSFORM  FOR  THE  1ST 
TRANSFORM  HALF  THE  DATA  AT 
HALF  BY  CONJUGATE  SYMMETRY. 

REAL  TRANSFORM  FOR  THE  1ST  DIMENSION,  N  EVEN.  METHOD — 
TRANSFORM  A  COMPLEX  ARRAY  OF  LENGTH  N/2  WHOSE  REAL  PARTS 
ARE  THE  EVEN  NUMBERED  REAL  VALUES  AND  WHOSE  IMAGINARY  PARTS 
ARE  THE  ODD  NUMBERED  REAL  VALUES.  SEPARATE  AND  SUPPLY 
THE  SECOND  HALF  BY  CONJUGATE  SYMMETRY. 
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OATA(  I3)=DATA{ J3) 

PFTT172 

DATA{ I3+1)=DATA( J3+1) 

FFTT173 

DATA( J3)=TEMPR 

FFTT174 

125 

DATA( J3+1)=TEMPI 

FFTT175 

130 

M=NP2HF 

FFTT176: 

140 

IF(J-M)150,150,145 

FFTT177 

145 

J  =  J-M 

FFTT17S. 

M=M/2 

FFTT179( 

IF(M-N0N2)150,140,140 

FFTT180l 

150 

J=J+M 

FFTT181 

C 

FFTT182' 

C 

MAIN  LOOP  FOR  FACTORS  OF  TWO.  PERFORM 

FOURIER  TRANSFORMS  OF  FFTT183 

C 

LENGTH  FOUR,  WITH  ONE  OF  LENGTH  TWO  IF 

NEEDED.  THE  TW 

IDDLE  FACT0RFFTT184 

C 

W=EXP{ ISIGN*2*PI*SQRT(-1)*M/(4«MMAX) ). 

CHECK  FOR  W=IS 

IGN’f'SQRT(-l)  FFTT185 

C 

AND  REPEAT  FOR  W= I S I GN*SQRT (-1 ) *CON J UG ATE ( W ) . 

FFTl 186, 

C 

FFTT1871 

N0N2T=N0N2+N0N2 

FFTT188; 

IPAR=NTW0/NP1 

FFTT189 

310 

IF( IPAR-2)350,330,320 

FFTT190 

320 

IPAR=IPAR/4 

FFTT191i 

GO  TO  310 

FFTT192 . 

330 

DO  340  I1=1,I1RNG,2 

FFTT193C 

DO  340  J3=I1,N0N2,NP1 

FFTT194i 

DO  340  K1=J3,NTOT,NON2T 

FFTT195V 

K2=K1+N0N2 

FFTT196, 

TEMPR=DATA{K2) 

FFTT197  . 

TEMPI=DATA(K2+1) 

FFTT198! 

DATA{K2) =DATA(K1)-TEMPR 

FFTT199! 

DATA! K2+1)=DATA(K1+1)-TEMPI 

FFTT200 

DATA! Kl) =DATAIK1)+TEMPR 

FFTT201> 

340 

DATA( K1+1)=DATA(K1+1)+TEMPI 

FFTT202. 

350 

MMAX=NCN2 

FFTT203. 

360 

I  F(MMAX-NP2HF >370,600, 600 

FFTT204 

370 

LMAX=MAX0{N0N2T,MMAX/2) 

FFTT205 

I F( MMAX-N0N2 >405,405, 380 

FFTT206 

380 

THET  A=-T  WOP  I  AFLOAT  (N0N2)/ FLOAT  (4=l=MMAX) 

FFTT207 

IFI ISIGN)400,390,390 

FFTT208 

390 

THETA=-THETA 

FFTT209; 

400 

WR=COS(THETA) 

FFTT2101 

WI=SIN(THETA) 

FFTT211 

WSTPR=-2.*WI*WI 

FFTT212( 

WSTPI  =2.*WR=I=WI 

FFTT213C 

405 

DO  570  L=N0N2,LMAX,N0N2T 

FFTT214C 

M=L 

FFTT215i 

I  F{MMAX-N0N2) 42  0,42  0,410 

FFTT216C 

410 

W2R=WR=!=WR-WI=4'WI 

FFTT217V 

W2I=2 .*WR*WI 

FFTT218C 

W3R=W2P-WR-W2I-WI 

FFTT219.: 
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W3I  =  W2R=«'WI+W2I*WR 
420  DO  530  Ii=l,IlRNG,2 

DO  530  J3=I1»N0N2,NP1 
K'1IM=J3+IPAR*M 
IF(MMAX-NON2)430,43C,440 
430  KMIN=J3 
440  KDIF= IPAR*MMAX 
450  KSTEP=4’!'KDIF 

DO  520  K1=KMIN,NT0T,KSTEP 

K2=K1+KDIF 

K3=K2+KDIF 

K4=K3+KDIF 

IF(MMAX-N0N2) 460,460,480 
460  U1R=0ATA(K1)+DATA(K2) 

U1I=DATA(K1+1)+DATA(K2+1) 
U2R=DATA(K3)+0ATA{K4) 
U2I=DATA(K3+1)+DATA(K4+1) 
U3R=0ATA(K1)-DATA(K2) 
U3I=DATA{K1+1)-DATA(K2+1) 
IF(ISIGN)470,475,475 
470  U4R=DATA(K3+1)-DATA(K4+1) 

U4I=DATA (K4)-DATA(K3) 

GO  TO  510 

475  U4R=DATA(K4+1)-DATA(K3+1) 

U4I=DATA(K3)-DATA(K4) 

GO  TO  510 

480  T2R=W2P*DATA{K2)-W2I=!'DATA(  K2+1 ) 
T2I=W2R*DATA(K2+1)+W2I*DATA(K2) 
T3R=WR*DATA{  K3  J-WI=!'DATA(K3  +  1) 
T3I=WR4DATA( K3+1 )+WI«DATA(K3} 
T4R=W3R*DATA(K4)-W3I*DATA(K4+1) 
T4I=W3R4DATA(K4+1 )+W3I*DATA(K4) 
U1R=DATA(K1)+T2R 
U1I=DATA(K1+1)+T2I 
U2R=T3R+T4R 
U2I=T3I+T4I 
U3R=DATA(K1)-T2R 
U3I=DATA(K1+1)-T2I 
IF(ISIGN)490,500,500 
490  U4R=T3I-T4I 
U4I=T4P-T3R 
GO  TO  510 
500  U4R=T4I-T3I 
U4I=T3R-T4R 

510  DATA( K1)=U1P+U2R 

DATA( K1+1)=U1I+U2I 
DATA( K2) =U3R+U4R 
DATA( K2+1) =U3I+U4I 
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FFTT223: 

FFTT221. 

FFTT222C 

FFTT223L 

FFTT224r 

FFTT225C 

FFTT226C 

FFTT227C 

FFTT228': 

FFTT229C 

FFTT230C 

FFTT231: 

FFTT232C 

FFTT233C 

FFTT234(, 

FFTT235C 

FFTT236: 

FFTT237C 

FFTT238; 

FFTT239C 

FFTT240C 

FFTT241C 

FFTT242C 

FFTT243C 

FFTT244C 

FFTT245C 

FFTT246C 

FFTT247C 

FFTT248C 

FFTT249C 

FFTT250C 

FFTT251C, 

FFTT252C 

FFTT253C 

FFTT254C 

FFTT255C 

FFTT256C 

FFTT257C 

FFTT258C 

FFTT259C 

FFTT260C 

FFTT26,1C 

FFTT262C 

FFTT263C 

FFTT264C 

FFTTZcSO 

FFTT266C 

FFTT267C 
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DATA( K3)=U1R-U2R 

FFTT263 

DATA(K3+1)=U1I-U2I 

FFTT269 

DATA{ K4)=U3R-U4R 

FFTT270 

520 

DATA{ K4+1)=U3I-U4I 

FFTT271 

KMIN=4*( KMIN-J3)+J3 

FFTT272. 

KDIF=KSTEP 

FFTT273> 

IF(KDIF-NP2)450,530,530 

FFTT274. 

530 

CONTINUE 

FFTT275( 

M=MMAX-M 

FFTT276( 

IF(ISIGN)540,550,550 

FFTT277 

540 

TEMPR=WR 

FFTT278i 

WR=-WI 

FFTT279(. 

WI=-TEMPR 

FFTT280 

GO  TO  560 

FFTT281i 

550 

TEMPR=WR 

FFTT282 . 

WR=WI 

FFTT283. 

W I  =T  E  M  PR 

FFTT284( 

560 

IF(M-LMAX)565,565,410 

FFTT285:; 

565 

TEMPR=WR 

FFTT286( 

WR=WR*WSTPR-WI4WSTPI+WR 

FFTT287. 

570 

WI=WI  ’i'WSTPR+TEMPR^WSTPI+WI 

FFTT288 

IPAR=3-IPAR 

FFTT289: 

MMAX=MMAX+MMAX 

FFTT290(: 

GO  TO  360 

FFTT291e 

C 

FFTT292C 

C 

MAIN  LOOP  FOR  FACTORS  NOT  EQUAL  TO  TWO.  APPLY 

THE  TWIDDLE  FACTOR 

FFTT293: 

c 

W  =  EXP(  ISIGN*2=i'PI*SQPT(-l)*(J2-l)*(Jl-J2)/{NP2*IFPl)  )  ,  THEN 

FFTT2941 

c 

PERFORM  A  FOURIER  TRANSFORM  OF  LENGTH  IFACT(IF) 

,  MAKING  USE  OF 

FFTT2951 

c 

CONJUGATE  SYMMETRIES. 

FFTT296 

c 

FFTT297i 

600 

I  F(NTW0-NP2)  605 ,700,700 

FFTT2981 

605 

IFP1=N0N2 

FFTT299 

IF  =  1 

FFTT300( 

NPlHF=NPl/2 

FFTT30L 

610 

IFP2= IFP1/IFACT( IF) 

FFTT302> 

J  1RNG  =  NP2 

FFTT303 

IF( ICASE-3)612, 611,612 

FFTT304: 

611 

J1RNG=(NP2+IFP1 )/2 

FFTT3051 

J2STP=NP2/IFACT(IF) 

FFTT306i 

JlRG2=(J2STP+IFP2)/2 

FFTT307:, 

612 

J2MIN=i+IFP2 

FFTT3031 

IF( IFP1-NP2) 615,640,640 

FFTT3091 

615 

DO  635  J2=J2MIN,IFP1,IFP2 

FFTT310V 

THETA=-TWOPI*FLOAT( J2-1)/ FLOAT (NP2) 

FFTT3il( 

IF{ ISIGN)625,620,620 

FFTT312 : 

620 

thfta=-theta 

FFTT313C 

625 

SINTH=SIN(THETA/2. ) 

FFTT3141 

WSTPR=-2.*SINTH*SINTH 

FFTT315 
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WSTPI=SIN(THETA) 

WR=WSTPR+1, 

WI=WSTPI 

J1MIN=J2+IFP1 

DO  635  J1  =  J1MIN, JlRNGt  IFPl 

IlMAX=Jl+IlRNG-2 

DO  630  I1=J1, IlMAXt2 

DO  630  13=11, NT0T,NP2 

J3MAX=I3+IFP2-NP1 

DO  630  J3=I3, J3MAX,NP1 

TEMPR=DATA( J3) 

DATA( J3) =DATA{ J3 ) *WR-OAT  A ( J3+1 ) -W I 
630  DATA(  J3+l)=TEMPR’i'WI+0ATA(  J3+1)*WR 
TEMPR=WR 

WR=WP*WSTPR-WI*WSTPI+WR 
635  WI=TEMPR=^WSTPI+WI=i=WSTPR+WI 
640  THETA=-TWOP I/FLOAT ( IF  ACT ( IF) ) 
IF{ISIGN)650,645,645 
645  THETA=-THETA 
650  SINTH=SIN(THETA/2. ) 

WSTPR=-2.*SINTH*SINTH 

WSTPI=SIN{THETA) 

KSTEP=2*N/IFACT( IF) 

KRANG=KSTEP#( IFACT( IF)/2)+l 

DO  698  I1=1,I1RNG,2 

DO  698  13=11, NT0T,NP2 

DO  690  KMIN=1,KRANG,KSTEP 

J1MAX=I3+J1PNG-IFP1 

DO  680  J1=I3 , JIMAX, IFPl 

J3MAX=J1+IFP2-NP1 

DO  680  J3=J1,J3MAX,NP1 

J2MAX=J3+IFP1-IFP2 

K=KMIN+( J3-J1+(J1-I3)/IFACT( IF))/NP1HF 
IF(KMIN-1)655,655,665 
655  SUMR=0. 

SUMI=0. 

DO  660  J2=J3 , J2MAX, IFP2 
SUMR=SUMR+DATA( J2) 

660  SUMI=SUMI+DATA( J2+1) 

WORK{K)=SUMR 
W0RK(K+1)=SUMI 
GO  TO  680 

665  KC0NJ=K+2*(N-KMIN+1 ) 

J2=J2MAX 
SUMR=DATA( J2) 

SUMI=DATA( J2+1 ) 

0L0SR=0. 

0LDSI=0. 


FFTT3160 

FFTT3170 

FFTT3180 

FFTT3190 

FFTT320D 

FFTT3210 

FFTT3220 

FFTT323'0 

FFTT3240 

FFTT3250 

FFTT3260 

FFTT3270 

FFTT3280 

FFTT3290 

FFTT3300 

FFTT3310 

FFTT3320 

FFTT3330 

FFTT3340 

FFTT3350 

FFTT3360 

FFTT3370 

FFTT3380 

FFTT3390 

FFTT3400 

FFTT3410 

FFTT3420 

FFTT3430 

FFTT3440 

FFTT3450 

FFTT3460 

FFTT3470 

FFTT3480 

FFTT3490 

FFTT350C 

FFTT3510 

FFTT352C 

FFTT3530 

FFTT3540 

FFTT3550 

FFTT356C 

FFTT3570 

FPTT3580 

FFTT3590 

FFTT3600 

FFTT3t!.  0 

FFTT3620 

FFTT3630 
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J2=J2-IFP2 
670  TEMPR=SUMR 
TEMPI=SUMI 

SUMR=TWOWR*SUMR-OLDSR+DATA( J2) 

SUMI=TWOWR*SUMI-OLOSI+DATA( J2+1) 

OLDSP=TEMPR 

OLDSI=TEMPI 

J2=J2-IFP2 

IF( J2-J3)675,675,670 
675  TEMPR=WR«SUMR-OLDSR+DATA( J2) 
TEMPI=WI*SUMI 
WORK! K)=TEMPR-TEMPI 
WORK( KCONJ)=TEMPR+TEMPI 
TEMPR  =  WR’!'SUMI-OLDSI+DATA(  J2  +  1) 
TEMPI=WI«SUMR 
WORK( K+1)=TEMPR+TEMPI 
WORK( KCONJ+1 )=TEMPR-TEMPI 
680  CONTINUE 

IF(KMIN-1)685,685,686 

685  WR=WSTPR+1. 

WI=WSTPI 

GO  TO  690 

686  TEMPR=WR 
WRsWR^WSTPR-WI’l'WSTPI+WR 
WI=TEMPR=i=WSTPI  +  WI*WSTPR+WI 

690  TWOWR=WR+WR 

IF( ICASE-3)692,691,692 

691  IF{ IFP1-NP2)695,692,692 

692  K=1 
I2MAX=I3+NP2-NP1 

DO  693  I2=I3tI2MAX,NPl 
DATA( I2)=W0RK(K) 

DATA! I2+1)=W0RK(K+1) 

693  K=K+2 

GO  TO  698 
C 

C  COMPLETE  A  REAL  TRANSFORM  IN  THE 

C  JUGATE  SYMMETRIES  AT  EACH  STAGE. 

C 

695  J3MAX=I3+IFP2-NP1 

DO  697  J3=I3, J3MAXtNP1 

J2MAX=J3+NP2-J2STP 

DO  697  J2=J3,J2MAX, J2STP 

J1MAX=J2+J1RG2-IFP2 

J1CNJ=J3+J2MAX+J2STP-J2 

DO  697  J1=J2,J1MAX,  IFP2 

K=1+J1-I3 

DATA( J1)=W0RK(K) 
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FFTT364 

FFTT365 

FFTT366 

FFTT367 

FFTT368V 

FFTT369 

FFTT3701 

FFTT371C 

FFTT372 ; 

FFTT373i 

FFTT374C 

FFTT375 

FFTT3  76'. 

FFTT377. 

FFTT;;78i. 

FFTT379( 

FFTT380 . 

FFTT381{ 

FFTT382< 

FFTT383 

FFTT3841 

FFTT385( 

FFTT386: 

FFTT3871 

FFTT3881 

FFTT389: 

FFTT390t 

FFTT3911. 

FFTT392: 

FFTT393V 

FFTT394 . 

FFTT395C 

FFTT396i 

FFTT397 . 

FFTT398’. 

FFTT399'; 

1ST  DIMENSION,  N  ODD,  BY  CON-  FFTT400l 

FFTT401" 
FFTT402( 
FFTT403 
•  FFTT404C 
FFTT405( 
FFTT406' 
FFTT407C 
FFTT408 ■ 
FFTT^OOC 
FFTT^-iO. 
FFTT411(; 
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696 

697 

698 


C 

C 

c 

c 

700 

701 


702 

703 


710 


720 


725 

730 

731 


2.0  FQURT  DAT 

DATA( Jl+1) =W0RK(K+1) 

IF( Jl-J2)697t697,696 
DATA( J1CNJ)=W0RK(K) 

DATA{ JlCNJ+1 )=-WORK(K+l ) 

J1CNJ=J1CNJ-IFP2 

CONTINUE 

IF=IF+1 

IFP1=IFP2 

IF( IFP1-NP1)700»700,610 

COMPLETE  A  REAL  TRANSFORM  IN  THE  1ST 
JUGATE  SYMMETRIES. 

GO  TO  {900,800,900»701) , ICASE 

NHALF=N 

N=N+N 

THETA=-TWOPI/FLOAT(N) 

IF{ ISIGN)703t702,702 
THETA=-THETA 
SINTH=SIN(THETA/2.  ) 
WSTPR=-2.«SINTH*SINTH 
WSTPI=SINITHETA) 

WR=WSTPR+1. 

WI=WSTPI 

IMIN=3 

JMIN=2*NHALF-1 
GO  TO  725 
J=JMIN 

DO  720  I=IMIN,NT0T,NP2 
SUMR=(DATA(I)+DATA(J) )/2. 

SUMI=( DATA! I+l )+DATA( J+1) )/2. 
DIFR=(CATA( I)-DATA( J) )/2. 

DIFI  =  ( DATA! I  +  l ) -DAT  A ( J  + 1 ) ) / 2 . 

TEMPR=WR*SUMI+W I*DIFR 

TEMPI =WI*SUMI-WR-DI FR 

DATA(  I)=SUMR+TEMPR 

DATA! I+1)=DIFI+TEMPI 

DATA! J)=SUMR-TEMPR 

DATA(J+1)=-DIFI+TEMPI 

J=J+NP2 

IMIN=IMIN+2 

JMIN= JMI N-2 

TEMPR=WR 

WR=WR=l=WSTPR-WI*WSTPI+WR 
WI=TEMPR*WSTPI+WI*WSTPR+WI 
IF( IMIN-JMIN) 710,730t740 
IF(ISIGN)731»740,740 
DO  735  I=IMI N,NT0T,NP2 
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735  DATA( I+l )=-DATA( I+l) 

740  NP2=NP2+NP2 

NTOT=NTOT+NTOT 
J=NT0T+1 
IMAX=NT0T/2+l 
745  IMIN=IMAX-2*NHALF 
I=IMIN 
GO  TO  755 

750  DATAIJ)=DATA(I  ) 

DATA( J+1 )=-DATA( I  +  l ) 

755  1=1+2 

J=J-2 

IF( I-IMAX)750,760,760 
760  DATA! J )=DATA( IM I N ) -DATA ( I M I N+1 ) 

DATA( J+1)=0. 

IF(I-J )770,780,780 
765  DATA( J)=OATA(I) 

DATA(  J  +  U  =  DATA(  I  +  l) 

770  1=1-2 

J=J-2 

IF( I-IMIN)775t775,765 
775  DATA! J)=DATA{ IMIN)+DATA( IMIN+1) 

DATA! J+1 )=0. 

IMAX=IMIN 
GO  TO  745 

780  ,DATA( 1)=DATA{ 1)+DATA(2) 

DATA! 2)=0. 

GO  TO  900 
C 

C  COMPLETE  A  REAL  TRANSFORM  FOP  THE  2N0  OR  3PD  DIMENSION  BY 

C  CONJUGATE  SYMMETRIES. 

C 

800  IF( IlRNG-NPl )805,900,900 
805  DO  860  I3=1,NT0T,NP2 
I2MAX=I3+NP2-NP1 
DO  860  12=13, I2MAX,NP1 
IMIN=I2+I1RNG 
IMAX=I2+NPl-2 
JMAX=2*I3+NP1-IMIN 
IF(I2-I3)820,820,810 
810  JMAX=JMAX+NP2 
820  IF( IDIM-2)850,850,830 
830  J=JMAX+NPO 

DO  840  I=IMIN, IMAX,2 
OATA( I )=DATA( J) 

DATA! I+l )=-DATA( J+1) 

840  J=J-2 

850  J=JMAX 


FFTT460C 

FFTT461C 

FFTT462C 

FFTT463C 

FFTT4640 

FFTT465C 

FFTT466C 

FFTT4670 

FFTT468C 

FFTT469C 

FFTT470C 

FFTT471C 

FFTT472C 

FFTT473C 

FFTT474C 

FFTT4750 

FFTT476C 

FFTT477C 

FFTT4730 

FFTT479C 

FFTT480C 

FFTT481C 

FFTT482C 

FFTT4830 

FFTT484C 

FFTT485G 

FFTT486G 

FFTT487C 

FFTT488C 

FFTT489: 

FFTT490C 

FFTT491C 

FFTT492C 

FFTT493C 

FFTT494.:. 

FFTT495t. 

FFTT496C 

FFTT497C 

FFTT498C 

FFTT499C 

FFTT500t 

FFTT501C 

FFTT502: 

FFTT503C 

FFTT504C 

FFTT505C 

FFTT:j06C 

FFTT507C 


o  o  o 


-  103  - 


RELEASE  2.0 


FOUR! 


DATE  =  77139  21/11/39 


DO  860  I=IMIN,IMAX,NPO 
DATA! I )=DATA( J) 

DATAi I+li=-DATA{ J+1) 

860  J=J-NPO 

END  OF  LOOP  ON  EACH  DIMENSION 

900  NP0=NP1 
NP1=NP2 
910  NPREV=N 
920  RETURN 
END 
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